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Executive  Summary 


The  research  described  in  this  report  was  conducted  in  support  of  Strategic  Environmental  Research  and  De¬ 
velopment  Program  (SERDP)  SEED  Broad  Agency  Announcement  (BAA),  Statement  of  Need  UXSEED-05-02, 
“Innovative  Approaches  to  Unexploded  Ordnance  (UXO)  Cleanup.”  A  SERDP  SEED  research  and  development 
project  UX-1446  entitled  “A  Unified  Approach  to  UXO  Discrimination  Using  the  Method  of  Auxiliary  Sources” 
was  proposed  in  response  to  the  above  BAA. 

The  main  emphasis  of  this  research  was  to  explore  the  fundamental  characteristics  of  the  Surface  Magnetic 
Charge  (SMC)  and  Standard  Excitation  Approach  (SEA)  methods  when  applied  to  UXO  discrimination  problems. 
Both  methods  were  derived  from  the  Method  of  Auxiliary  Sources  (MAS),  and  thus  represent  the  secondary 
magnetic  fields  from  a  compact  metallic  target  with  a  surface  of  magnetic  charge.  The  SEA  and  SMC  are  relatively 
new  modeling  techniques  for  UXO  discrimination.  Therefore,  we  investigated  some  fundamental,  as  well  as 
practical,  characteristics  of  the  forward  model.  These  include  the  accuracy  with  which  the  methods  can  model 
sensor  data,  the  speed  to  carry  out  the  forward  modeling,  and  the  type  of  discrimination  algorithms  amenable  to 
each  of  the  forward  modeling  methods.  For  the  SEA,  we  wanted  to  determine  the  ease  with  which  the  sources  can 
be  derived  for  a  particular  target.  For  the  SMC,  we  wanted  to  determine  if  the  surface  magnetic  charge  distribution 
is  a  good  discriminant,  and,  if  so,  what  algorithm  is  required  to  obtain  a  stable  estimate  of  the  magnetic  charge. 

Several  data  sets  were  analyzed  as  part  of  this  investigation.  Geonics  EM63  data  were  collected  on  the  Engi¬ 
neer  Research  and  Development  Center  (ERDC)  test  stand  in  Vicksburg,  MS  by  Sky  Research  Inc.,  University  of 
British  Columbia  personnel,  and  Cliff  Morgan  and  Morgan  Field  of  the  U.S.  Army  Corps  of  Engineers  (USACE)- 
ERDC  from  March  16-April  11,  2006.  Geophex  GEM3  data  were  collected  on  the  ERDC  test  stand  in  Vicksburg 
MS  by  Sky  Research  Inc.,  and  Cliff  Morgan  and  Morgan  Field  of  the  USACE-ERDC  from  February  2-10,  2006. 
Geophex  GEM3  data  and  Geonics  EM63  data  collected  during  the  month  of  July,  2005  on  the  Sky  Research  Inc. 
UXO  testplot  in  Ashland,  OR  were  also  processed  for  this  project.  In  this  report,  we  do  not  estimate  location 
and  orientation  from  the  data.  Instead,  we  assume  that  the  location  and  orientation  estimates  can  be  provided  by 
alternative  means. 
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Chapter  1 

Introduction 


The  research  described  in  this  report  was  conducted  in  support  of  SERDP  SEED  Broad  Agency  Announcement, 
Statement  of  Need  UXSEED-05-02,  “Innovative  Approaches  to  UXO  Cleanup.”  A  SERDP  SEED  research  and 
development  project  UX-1446  entitled  “A  Unified  Approach  to  UXO  Discrimination  Using  the  Method  of  Auxil¬ 
iary  Sources”  was  proposed  in  response  to  the  above  BAA. 

Several  data  processing  techniques  for  geophysical  survey  data  have  been  developed  for  discriminating  be¬ 
tween  UXO  and  non-UXO  items.  Current  Electromagnetic  Induction  (EMI)  data  processing  practices  involve 
inverting  sensor  data  for  dipole  parameters,  and  then  applying  statistical  classification  techniques  to  the  recovered 
parameters.  The  dipole  model  is  an  attractive  modeling  technique  for  UXO  discrimination  because:  (1)  it  is  very 
fast,  (2)  sensor  data  collected  over  UXO  are  nearly  dipolar  if  we  are  not  measuring  the  field  at  a  close  distance  to 
the  target  (relative  to  target  size),  and  (3)  the  parameters  of  the  dipole  model  have  a  good  physical  basis  and  we 
have  an  understanding  of  how  material  properties,  size,  and  shape  affect  the  dipole  model  parameters.  However, 
the  secondary  field  for  complex,  heterogeneous  targets  will  often  appear  non-dipolar,  particularly  if  we  are  close 
to  the  target.  Figure  1.1  compares  the  time  domain  electromagnetic  response  for  a  horizontal  105  mm  projectile 
that  is  aligned  parallel  to  the  y-axis.  The  data  were  measured  1  m  above  the  target,  using  a  Geonics  EM63  time 
domain  electromagnetic  sensor.  Even  though  the  data  are  collected  1  m  from  the  target,  there  are  still  non-dipolar 
components  in  the  signal.  If  the  response  was  truly  dipolar,  we  would  expect  a  dipole  whose  axial  excitation  was 
horizontal  and  the  response  would  have  symmetry  with  respect  to  both  x  and  y  axes. 

Typical  EMI  data  processing  steps  are  to  invert  sensor  data  for  dipole  parameters,  and  then  to  apply  statistical 
classification  techniques  to  the  recovered  parameters.  A  dipole  parameter  inversion  method  minimizes  a  data 
misfit  function,  and  therefore  attempts  to  fit  both  the  dipolar  and  non-dipolar  components  of  the  signal.  In  the 
context  of  a  dipole  inversion,  the  non-dipolar  components  of  the  secondary  field  are  essentially  data  noise.  As 
a  result  this  data  noise  is  propagated  through  the  inversion  to  produce  inaccuracies  in  the  recovered  parameters. 
The  model  parameters  (i.e.,  the  position,  orientation  and  polarization  tensor  components)  will  adjust  to  fit  the 
non-dipolar  components  of  the  data.  For  example,  the  dipole  inversion  applied  to  Figure  1 . 1  predicts  a  target  with 
a  5.7  degree  dip  in  order  to  compensate  for  the  non-symmetry  in  the  data. 

The  initial  objective  of  this  research  project  was  to  apply  the  MAS  as  a  modeling  technique  for  inversion  of 
UXO  sensor  data.  The  inputs,  or  model  vector,  for  the  MAS  would  be  the  target  size,  shape,  and  electromagnetic 
material  properties.  The  MAS  can  model  the  non-dipolar  components  of  the  secondary  field.  This  suggests  that  the 
spread  of  feature  parameters  for  a  MAS  model  would  be  relatively  smaller  than  for  dipole  parameters.  The  overall 
objective  of  the  project  was  to  develop  a  general  inversion  methodology  based  on  the  MAS  that  could  be  used  to 
invert  any  combination  and  configuration  of  magnetic.  Transient  Electromagnetic  Method  (time  domain)  (TEM) 
and  Finite  Element  Method  (frequency  domain)  (FEM)  induction  data.  Specific  objectives  included  (1)  determine 
the  feasibility  of  using  the  MAS  forward  model  as  the  basis  for  the  inversion  of  total-field  magnetics,  TEM  and 
FEM  induction  data;  (2)  investigate  the  feasibility  of  using  the  MAS  forward  model  as  the  basis  for  joint  and/or 
cooperative  inversion  of  any  combination  of  magnetic,  TEM  and  FEM  induction  data;  and  (3)  determine  the 
discrimination  capability  of  the  physical  parameters  recovered  by  our  inversion  for  different  combinations  and 
configurations  of  sensors. 

Through  numerous  comparisons  between  measured  data,  analytical  solutions,  and  data  modeled  by  the  MAS, 
we  found  that  the  MAS  is  able  to  accurately  model  compact  metallic  objects,  such  as  UXO.  The  MAS  accurately 
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(a)  Image  of  the  20th  time  channel  of  EM63  data 


(b)  Profile  directly  over  the  center  of  the  target  (x  =  3.59  m) 

Multiple  channels  for  line  10  Channel  ’residual’  for  line  10 


(c)  Recovered  Parameters 


X  =  3.60  (±0.00)  m 
Y  =  2.59  (±0.03)  m 
Z  =  -0.96  (  ±0.22)  m 
0  =  179.9°  (±0.5°)  deg 
0  =  94.3°  (±0.2°)  deg 
CorrCoeff  =  0.975 


Figure  1.1:  Results  from  inverting  Geonics  EM63  data  collected  over  an  105  mm  projectile. 


modeled  the  EMI  response  over  a  large  frequency  range  that  extended  to  the  static  case  (Figure  1 .2,  from  Beran, 
2005,  Chapter  2).  The  ability  to  model  over  such  a  large  range  made  the  MAS  an  attractive  modeling  technique 
for  joint  inversion  of  multiple  data  sets.  However,  the  computational  time  required  by  the  MAS  to  calculate  sensor 
data  prevented  us  from  applying  the  MAS  for  inversion  purposes.  Considering  the  accuracy  of  the  MAS  method, 
we  revised  the  objective  of  this  study  to  evaluate  the  potential  of  two  modeling  techniques  derived  from  the  MAS 
approach:  the  Standardized  Excitation  Approach  and  the  Surface  Magnetic  Charge  approach.  Although  the  SEA 
and  SMC  are  both  surface  magnetic  charge  distribution  techniques,  they  represent  two  fundamentally  different 
approaches  to  UXO  discrimination.  The  SEA  is  used  as  part  of  a  library  or  hypothesis  testing  technique,  and  the 
SMC  is  part  of  a  parameter  estimation/statistical  classification  technique. 

The  SEA  approach  involves  determining  a  set  of  sources  that  are  a  function  of  target  properties  (i.e.,  shape, 
size,  conductivity,  and  magnetic  permeability)  only.  These  sources  are  determined  numerically  through  the  ap¬ 
plication  of  MAS  combined  with  a  Thin  Skin  Approximation  (TSA).  The  sources  can  then  be  forward  modeled 
to  generate  either  frequency  or  time  domain  sensor  data.  The  SMC  models  the  field  with  a  magnetic  charge  dis¬ 
tribution  on  a  fictitious  surface  that  encloses  the  target.  The  forward  modeling  is  fast  and,  therefore,  sensor  data 
can  quickly  be  inverted  for  the  charge  distribution.  For  a  given  location  and  orientation,  the  determination  of  the 
charge  distribution  is  a  linear  inverse  problem.  The  sum  of  the  charge  distribution,  i.e.,  the  total  surface  magnetic 
charge,  can  be  used  for  discrimination. 

This  report  is  organized  in  the  following  manner.  Chapter  Two  outlines  our  research  and  development  of 
the  SEA  approach.  The  SEA  approach  involves  determining  a  set  of  sources  that  represent  target  parameters 
that  can  be  forward  modeled  to  generate  sensor  data.  These  sources  can  be  calculated  for  a  number  of  targets  in 
order  to  generate  a  source  library.  Discrimination  is  achieved  by  determining  which  target  within  the  library  has 
the  greatest  likelihood  of  producing  the  anomaly.  Chapter  Three  outlines  our  research  into  the  SMC  modeling 
technique.  The  SMC  approach  models  sensor  data  with  the  field  from  a  magnetic  charge  distribution  on  a  fictitious 
surface  that  encloses  the  target.  The  forward  modeling  is  very  fast  and,  therefore,  sensor  data  can  be  inverted  for 
the  charge  distribution.  The  sum  of  the  charge  distribution,  i.e.,  the  total  surface  magnetic  charge,  can  be  used 
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Figure  1.2:  Modeling  magnetostatic  data  with  the  Method  of  Auxiliary  Sources  method.  A  spheroid  is  modelled 
using  the  Earth’s  magnetic  field  at  Yuma,  Arizona.  The  analytical  solution  for  a  spheroid  is  used.  The  Method  of 
Auxiliary  Sources  solution  was  obtained  by  assuming  a  frequency  of  Hr  6  Hz. 


for  discrimination.  Additional  work  studying  the  iso-sphere  surfaces  and  a  single  horizontal  charge  layer  is  not 
presented  here,  as  we  chose  to  focus  on  the  properties  of  the  total  magnetic  charge  applied  to  spheroidal  surfaces. 
Chapter  Four  summarizes  this  report  and  suggests  future  research  directions. 
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Chapter  2 


Investigation  of  EMI  Sensing  for  Metallic 
Objects:  a  Unified  and  Versatile 
Standardized  Excitation  Approach  in 
Frequency-  and  Time-  domain 

2.1  Abstract 

This  chapter  investigates  the  SEA  combined  with  the  MAS/TSA.  This  approach  is  a  numerical  technique  for 
computing  the  EMI  response  from  a  three-dimensional,  electromagnetically  heterogeneous  object  in  both  near 
and  far  fields.  The  objective  of  the  SEA  is  to  determine  a  set  of  characteristic  sources,  called  the  Reduced  Set  of 
Sources  (RSS),  associated  for  each  UXO.  These  sources  can  then  be  used  for  fast  modeling  of  the  EMI  response. 
The  full  EMI  solution  is  obtained  by  the  superposition  of  responses  to  the  spheroidal  excitation  modes.  In  this 
investigation,  we  determine  the  effective  spheroidal  modes  by  decomposing  the  primary  field  simulated  for  a 
loop  transmitter.  It  is  shown  that  the  spheroidal  series  can  be  safely  truncated  at  the  maximum  mode  numbers 
(M,  N)  =  (3, 8)  for  the  sensor  target  geometries  encountered  in  the  UXO  problem.  For  this  study,  the  frequency 
domain  RSS  has  been  determined  for  a  collection  of  UXO  items.  This  library  of  sources  has  been  extended  to  time 
domain  data  through  a  convolution  algorithm  that  can  use  an  arbitrary  current  waveform  in  a  TEM  system.  The 
experiment-based  results,  which  are  in  both  real  frequency-domain  (GEM-3)  and  time-domain  (EM63)  system, 
confirm  its  applicability  to  real-world  sensors.  As  another  illustration  of  using  a  RSS-library,  we  have  conducted 
an  identification  test  of  UXO  based  on  test  stand  data.  In  this  test,  we  demonstrate  that  identifying  a  candidate 
UXO  can  be  achieved  through  a  straightforward  pattern  matching  procedure  by  identifying  a  best  misfit  value 
below  a  threshold.  We  have  also  begun  formulating  the  inverse  problem  of  determining  an  object  location  and 
orientation  by  deriving  the  analytical  expressions  of  the  sensitivity  coefficients  of  the  fields  with  respect  to  these 
parameters.  A  simple  example  of  location  and  orientation  inversion  is  presented  here. 

2.2  Introduction 

Improved  cleanup  of  UXO  requires  data  processing  algorithms  for  discriminating  objects  of  interest  from  clut¬ 
ter  and  metallic  scrap.  Reliable  discrimination  minimizes  excavations  and,  therefore,  reduces  the  cost  and  any 
associated  danger.  Accurate  modeling  of  EMI  responses  from  buried  targets  is  an  important  component  of  a  dis¬ 
crimination  algorithm.  Several  forward  modeling  techniques  have  been  proposed.  The  most  common  approach  is 
to  approximate  a  target’s  EMI  response  with  a  point  dipole  forward  model.  The  dipole  approximation  is  appropri¬ 
ate  for  modeling  fields  in  the  far-held  regime,  and  has  been  successfully  employed  in  UXO  discrimination  (Geng 
et  al.,  1999;  Gao  et  al.,  2000;  Carin  et  al.,  2001;  Miller  et  al.,  2001;  Zhang  et  al.,  2003;  Pasion  and  Oldenburg, 
2001).  However,  the  near  held  EMI  response  often  contains  non-dipolar  components.  In  particular,  the  EMI 
response  near  a  heterogeneous  UXO  (consisting  of  multiple,  different  sections  such  as  head,  body,  tail  and  hns. 
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and  rotation  copper  band)  depend  strongly  on  which  parts  of  the  object  are  closer  to  the  sensor  and  on  the  degree 
of  coupling  between  different  parts.  In  such  a  case  the  target  would  be  poorly  modeled  with  a  point  dipole.  The 
importance  of  interaction  effects  when  simulating  target  responses  was  investigated  in  Shubitidze  et  al.  (2003); 
Shubitidze,  O’Neill,  Shamatava,  Sun  and  Paulsen  (2004);  Shubitidze,  O’Neill,  Sun  and  Shamatava  (2004).  One 
technique  that  explicitly  models  a  target’s  shape,  volume,  EM  properties  in  modeling  procedures  of  representative 
techniques  is  the  MAS  (Shubitidze  et  al.,  2002;  Shubitidze,  O’Neill,  Sun  and  Paulsen,  2004),  later  improved  with 
the  incorporation  of  the  TSA  (Sun  et  al.,  2002,  2004).  The  hybrid  MAS/TSA  technique  (Shubitidze,  O’Neill, 
Sun,  Shamatava  and  Paulsen,  2004)  does  not  need  fine  mesh  discretization  within  a  target  or  on  its  surface,  as 
stressed  in  the  other  numerical  schemes  like  the  Surface  Impedance  Boundary  Condition  (SIBC)-based  three  di¬ 
mensional  Methods  of  Moments  (MoM)  (Sebak  et  al.,  1991)  and  the  finite  difference  and  finite  element  methods 
(FDM/FEM)  when  dealing  with  the  situation  of  very  small  skin  depth.  These  mesh-dependent  techniques  can  be 
slow  and  even  inaccurate,  thus  making  them  impractical  in  this  context.  In  contrast,  the  MAS/TSA  technique  has 
demonstrated  an  ability  to  accurately  solve  the  Three-dimensional  (3D)  full  EMI  response  for  a  highly  conducting 
and/or  ferrous  object  in  a  relatively  short  time. 

A  class  of  SEA  (Braunisch  et  al.,  2001;  Ao  et  al.,  2002;  Barrows  et  al.,  2004;  Chen  et  al.,  2004;  Sun  et  al., 
2005;  Shubitidze  et  al.,  2005b)  has  been  recently  proposed  to  synthesize  the  full  EMI  solutions  by  the  superpo¬ 
sition  of  responses  to  the  fundamental  excitation  modes.  Importantly,  these  responses  are  unique  for  each  UXO 
and  can  be  stored  as  a  library  that  is  used  to  rapidly  produce  the  complete  EMI  solutions  for  any  location  and 
orientation  of  the  target.  Chen  et  al.  (2004)  and  Sun  et  al.  (2005)  applied  the  SEA  to  UXO  discrimination  by  a 
pattern  matching  scheme  using  the  library  for  some  UXO.  In  their  approach,  an  inverse  problem  was  solved  in 
which  a  source  library  is  determined  from  high  quality  sensor  data  collected  over  each  candidate  UXO  at  different 
distances  and  orientations.  The  accuracy,  reliability,  and  frequency  band  of  the  established  library  is  limited  by 
the  measurement  noise  and  numerical  errors  arising  from  the  inherent  ill-conditioning  of  the  problem.  Shubitidze 
et  al.  (2005b)  presented  an  alternative  way  to  implement  the  SEA  based  on  the  hybrid  MAS/TSA  technique.  In 
this  scheme,  the  modal  response  coefficients  (or  RSS)  are  distributed  on  a  fictitious  spheroid  enclosing  an  object. 
The  RSS  is  determined  by  employing  a  physically  complete  numerical  simulation  of  an  object’s  response  to  each 
fundamental  mode.  Determining  the  RSS  using  simulations  allows  us  to  avoid  the  non-uniqueness  and  ill-posed 
nature  of  fitting  noisy  sensor  data.  In  principle,  the  RSS  can  be  obtained  very  accurately  for  any  excitation  mode. 

Our  study  of  the  SEA  involves  three  major  issues.  First,  it  is  known  that  the  spheroidal  approach  theoretically 
involves  an  infinite  series  of  fundamental  spheroid  modes  for  synthesizing  fields.  However  it  has  to  be  truncated 
in  practice.  The  determination  of  the  effective  spheroid  modes  in  fact  is  fundamental  for  achieving  computational 
efficiency  while  maintaining  accuracy.  We  study  this  issue  using  simulated  primary  fields  from  loop  transmitters 
with  geometries  typical  for  UXO  surveys.  Second,  we  illustrate  that  the  modal  response  coefficients  of  a  target 
can  be  used  to  accurately  model  both  frequency  and  time  domain  sensor  data  by  presenting  the  tests  against  the 
measurements  collected  by  a  GEM-3,  a  wideband  frequency-domain  sensor  developed  by  Geophex  Ltd.,  and 
EM63,  a  time-domain  instrument  developed  by  Geonics  Ltd.  These  tests  are  used  to  demonstrate  the  capability  of 
the  MAS/TSA-based  spheroidal  approach  as  a  unified  EMI  modeling  tool  in  both  frequency  and  time  domain.  A 
convolution  algorithm  is  used  to  model  the  time  domain  response  from  an  arbitrary  transmitter  waveform.  Third, 
we  have  started  to  build  a  RSS-library  for  several  UXO  utilizing  the  MAS/TSA-based  SEA.  To  demonstrate 
the  use  of  a  RSS-library  at  this  stage,  we  apply  a  straightforward  pattern  matching  procedure  to  carry  out  the 
identification  tests  of  UXO  by  inspecting  a  minimum  misfit  value  between  the  measured  and  calculated  patterns 
of  secondary  fields. 


2.3  The  Standardized  Excitation  Approach  Method 

This  section  briefly  describes  the  standardized  (spheroidal)  excitation  approach  detailed  in  Shubitidze  et  al. 
(2005b). 

2.3.1  Decomposing  a  Primary  Field 

Under  the  quasi-magnetostatic  approximation,  a  magnetic  field  outside  of  an  object  is  irrotational.  It  follows 
that  the  corresponding  scalar  magnetic  potential  satisfies  the  Laplace  equation  and  can  be  expanded  in  a  Fourier- 
Legendre  series  (Arfken,  1995)  in  a  prolate  spheroidal  coordinates.  The  prolate  spheroid  is  chosen  since  it  can 
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Figure  2.1:  The  prolate  spheroidal  coordinate  system  (77,  £,(/>)  with  major  and  minor  semi-axes  of  a  and  b,  the 
interfocal  distance  of  d  =  2 \J a 2  —  b2,  and  £  =  £0  =  2a /d. 


conform  to  an  elongated  body  of  revolution,  a  typical  geometry  of  UXO.  Let  us  denote  (77,  £,  <f>)  as  the  standard 
prolate  spheroidal  coordinates  with  —  1  <  77  <  1,  1  <  £  <  00,  0  <  <fi  <  2i r  and  the  inter-focal  distance  of  d,  as 
shown  in  Figure  2.1.  On  a  fictitious  spheroid  £  =  £0  surrounding  an  object,  a  primary  potential  field  Lpr  can  be 
expressed  as: 


,  1V1  IV  1 

rr(v,z 0,0)  =  -f-  E 

ra=0  n—m  p— 0 


(2.1) 


where  P™  are  associated  Legendre  functions  of  the  first  kind  with  degree  n  and  order  m  (Arfken,  1995),  and 
Tprn(4>)  is  cos{m4>)  forp  =  0  and  is  sin(m(f f>)  forp  =  1,  bpmn  the  spheroidal  expansion  coefficients.  Eq.  (2.1)  is  a 
decomposition  of  a  primary  magnetic  potential  into  a  number  of  fundamental  excitations  P™(p)P™(£)Tpm((/))  = 
qPprnn(rl-  £o;  <'/')  with  the  maximum  numbers  of  M  and  N,  being  infinity  in  theory.  Taking  a  gradient  operation 
of  Eq.  (2.1),  one  namely  obtains  the  decomposition  of  the  primary  magnetic  field  Hpr.  Its  normal  component 
fT|r(77,  £,  </>)  on  the  spheroid  can  be  written  as 


H0d  1 

~^TYi 


M  N  1 

m=0  n=m  p= 0 


(2.2) 


where  /p'm  is  the  first-order  derivative  of  the  associated  Legendre  function  with  respect  to  £,  is  the  metrical 
coefficient.  By  the  orthogonality  of  the  associated  Legendre  functions,  the  spheroidal  expansion  coefficients  bpmn 
can  be  derived  as 

2tj  4-  1  (n  —  mV  I"1  f27T 

bp mn  =  ~T THodP!r(lio)  (n  +  m)\ l  HfirlAoA)TPm{4>)dcPdr1,  (2.3) 

where  7  =  2  for  m  =  p  =  0  and  7=1  otherwise.  The  integration  in  (2.3)  is  evaluated  by  a  numerical  integration. 
This  completes  the  decomposition  of  a  primary  field  quantity  under  quasi-static  approximation. 
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Figure  2.2:  MAS/TSA-based  SEA.  Reduced  set  of  sources  gfmn  distributed  along  rings  on  the  spheroid  surface 
and  the  auxiliary  sources  Q\""‘  distributed  on  an  internal  auxiliary  surface  of  a  cylinder-like  object.  Both  induced 
sources  are  the  response  to  fundamental  excitation  fields  ippmn  and  produce  the  same  scattered  field  V'pmn  outside 
of  the  object. 

2.3.2  Synthesizing  the  Secondary  Fields 

We  are  interested  in  dealing  with  a  Body  of  Revolution  (BOR)  that  has  rotational  symmetry  about  the  z-axis  of  a 
Cartesian  coordinate  in  Figure  2.1.  The  symmetry  property  implies  that  the  primary  fields  behave  in  a  azimuthal 
dependence  of  Tpm((f>)  as  a  Fourier  series  (as  seen  in  Eqs.  (2.1)  and  (2.2),  and  each  primary  azimuthal  mode  in  the 
series  is  orthogonal.  As  a  result,  the  secondary  fields  also  preserve  this  angular  dependence  and  can  be  expressed 
similarly  in  a  Fourier  series  (Andreasen,  1965).  Let  us  consider  the  object’s  response  to  each  unit  excitation  H*”rm 
or  ifipmn-  We  can  introduce  the  equivalent  induced  magnetic  charges  distributed  along  the  ith  ring  on  the  fictitious 
(or  inside)  spheroidal  surface  in  light  of  the  rotational  symmetry  in  Fig.  2.2.  The  secondary  held  Hsc(r)  due  to 
an  object  can  be  written  as  a  linear  superposition  of  the  object’s  response  for  each  pmn  excitation  mode  (Chen 
et  al.,  2004;  Sun  et  ah,  2005),  i.e., 

M  N  1  Nq 

Hsc(r)  =  E  E  E  w  E  ?rnG(r,  r'),  (2.4) 

m—0  n=m  p—0  i=l 

where  r  is  the  position  vector  of  an  observation  point  outside  of  the  object  (the  spheroid),  the  r'  is  the  position 
vector  of  the  gfmn  along  the  ith  ring,  and  G(r.  r')  is  the  modal  Green’s  function  for  the  magnetic  held  and  is 
given  by 


G(r,  r') 


4-TTflo 


r  -  r: 


■'13 


PiTpm{(t>')d(j)' , 


(2.5) 


where  pi  is  the  radius  of  the  ith  ring.  The  modal  Green’s  function  can  be  understood  physically  as  the  vectorial 
held  radiated  by  a  unit  scalar  ring  source  harmonically  varying  in  strength  along  the  ring.  By  decomposing  the 
primary  helds  in  spheroidal  coordinates,  we  solve  the  scattering  problem  on  a  mode -by-mode  basis  and  then 
synthesize  the  scattered  held  by  applying  the  superposition  principle,  as  indicated  in  Eq.  (2.4). 

The  pattern  of  the  scattered  held  depends  on  intrinsic  attributes  of  the  scattering  object  (i.e.,  an  object’s 
geometry  and  physics),  as  well  as  extrinsic  attributes  of  the  object,  such  as  object  location,  orientation  and  an 
excitation  held.  One  way  to  look  at  UXO  discrimination  is  as  a  procedure  in  which  we  try  to  separate  the 
intrinsic  and  extrinsic  attributes.  From  Eqs.  (2.3)  and  (2.4),  extrinsic  attributes  are  characterized  by  the  spheroidal 
expansion  coefficients  bpmn  and  intrinsic  attributes  are  characterized  by  the  modal  response  amplitudes  qfmn. 
This  is  the  essence  of  the  SEA.  Therefore  once  such  response  amplitudes  are  found  for  an  object,  they  can  be 
stored  in  a  library  for  subsequent  rapid  computation  of  the  response  required  in  a  prospective  survey. 
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2.3.3  Solving  for  the  Modal  Response  Coefficients:  Reduced  Set  of  Sources 

Equation  (2.4)  suggests  that  a  straightforward  way  to  determine  the  RSS  is  to  form  an  inverse  problem  and  solve 
qPmn  eac|1  moc|e  corresponding  to  each  ring,  from  the  measured  data.  The  success  of  this  procedure  would 
require  numerous  independent  measurements  of  an  object  to  physically  reduce  the  degree  of  ill-posedness  which 
can  degrade  the  solution.  Recently,  Sun  et  al.  (2005)  applied  such  a  data-derived  approach  to  extracting  the 
modal  response  coefficients  for  each  candidate  by  carefully  designing  the  measurements  at  different  distances  and 
orientations.  However,  the  accuracy  and  reliability  of  the  model  parameters  determined  in  this  way  is  dependent 
on  the  measurement  noise  and  numerical  errors  arising  from  the  inherent  ill-conditioning  of  the  problem.  In 
addition,  a  library  generated  using  this  approach  is  limited  to  the  bandwidth  of  the  sensor. 

Shubitidze  et  al.  (2005b)  proposed  a  procedure  in  which  the  RSS  are  determined  by  employing  a  numerical 
simulation  of  an  object’s  response  to  each  fundamental  excitation  mode.  For  a  unit  excitation  ippmn  (Fig-  2.2),  the 
scattered  potential  field  outside  the  target  is  irrotational  and  can  be  written  as 

N. q 

Vpcmn(r)  =  ri),  (2-6) 

i—1 


where 


d(  r,r') 


1  1 
4tt  |  r  —  r'  |  ’ 


(2.7) 


is  the  potential  at  the  r  observation  point  produced  by  the  ith  ring  source  at  r'  on  the  spheroid  surface  where  the 
qPmn  are  located.  Eq.  (2.6)  assumes  that  the  field  V’pmn(r)  can  be  quite  accurately  expressed  for  observation 
points  outside  an  enclosing  fictitious  spheroid  using  q^mn.  For  a  given  candidate  target,  the  MAS  can  be  used 
to  compute  the  scattered  fields.  As  an  example,  let  us  consider  a  cylindrical  object  enclosed  by  a  spheroid  (Fig¬ 
ure  2.2).  The  MAS  distributes  a  set  of  auxiliary  charge  with  amplitudes  of  Qkmn,  k  =  1,  •  •  •  ,  Nq,  corresponding 
to  unit  excitation  mode,  on  an  auxiliary  surface  inside  the  scattering  object.  These  auxiliary  sources  are  obtained 
by  enforcing  the  EM  boundary  conditions  on  the  physical  surface  and  these  sources  are  responsible  for  the  exter¬ 
nal  fields  (Shubitidze  et  al.,  2002;  Shubitidze,  O’Neill,  Sun  and  Paulsen,  2004).  Thus  with  the  known  auxiliary 
charges  Qkmn,  the  scattered  potential  field  for  the  prriri  mode  can  be  written  as 


Nq 

V’pmnO)  =  X!  <9fem”flr(r>  rfe)’  (2-8) 

i=l 

where  v'k  is  the  position  vector  of  the  fcth  auxiliary  magnetic  charge  Qkmn,  located  inside  the  object.  The 
MAS  (Shubitidze  et  al.,  2002;  Shubitidze,  O’Neill,  Sun  and  Paulsen,  2004)  is  used  to  produce  the  synthetic 
observed  scattered  potentials  on  the  right-hand-side  of  Eq.  (2.6).  The  modal  response  coefficients  for  each  mode 
prrm  can  be  solved  via  a  linear  system  of  equations  Aq  =  %p  ,  where  A  is  J  x  Nq  a  matrix  whose  j ith  el¬ 
ements  given  by  Eq.  (2.7),  and  q  =  [q{mn,  q%mn,  ■  ■  ■  ,q1^nn]T,  a  TV^-dimensional  column  model  vector  and 

ip  =  [V’pmn(ri))  V’pmn(r2)i ' '  ■  ,  ippmn (r j) , ]T,  a  J-dimensional  column  data  vector.  Note  that  the  modal  re¬ 
sponse  coefficients  are  collectively  called  a  reduced  set  of  sources  because  its  number  is  substantially  less  than 
the  number  of  the  auxiliary  sources  Qkmn  required  in  the  MAS  (Shubitidze  et  al.,  2005b). 

This  model-based  approach  has  the  advantage  over  data-based  approaches,  because  the  corresponding  in¬ 
verse  problem  is  more  well-posed  and  the  RSS  can  be  more  accurately  obtained  for  any  excitation  mode  and  in 
an  ultra-wide  band  with  noise-free  synthetic  observed  data. 

2.4  Determination  of  Effective  Spheroidal  Modes 

As  described  in  section  2.3,  the  expressions  (2.1)  and  (2.2)  for  the  primary  and  secondary  fields  involve  infinitely 
many  terms  in  theory.  In  practice,  it  has  to  be  truncated  for  computational  efficiency.  Therefore,  to  implement 
the  SEA  it  is  important  to  find  the  suitable  maximum  mode  number  M  and  N  that  can  be  used  to  accurately 
represent  fields.  This  can  be  done  sufficiently  by  examining  the  decomposition  of  primary  fields;  modes  that 
have  a  negligible  contribution  to  the  primary  field  will  have  a  negligible  effect  on  the  secondary  fields.  Mode 
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Sensor 


Figure  2.3:  Sensing  of  a  BOR-like  object. 


truncation  requires  correct  determination  of  the  corresponding  spheroidal  expansion  coefficients  bpmn  for  each 
possible  mode.  In  this  investigation,  we  chose  to  decompose  uniform  and  non-uniform  fields  from  a  simulated 


GEM-3. 


2.4.1  Uniform  Primary  Field 

Let  us  consider  a  uniform  illuminating  field  H  =  z  and  decompose  the  field  on  a  spheroid  with  semi-major  and 
-minor  axes  of  a  =  11  cm  and  6  =  5  cm.  The  orientation  of  the  spheroid  is  described  by  the  azimuthal  and  polar 
angles  (a,  9),  as  shown  in  Fig.  2.3.  Due  to  rotational  symmetry,  we  can  illustrate  the  mode  expansion  to  represent 
H^1  along  a  generating  arc  defined  by  an  arbitrarily  chosen  spheroidal  azimuthal  coordinate  of  0  =  tt/4. 

Fig.  2.4  shows  the  comparison  between  the  exact  values  and  mode  expansions.  When  (M)  =  (o°,o°), 
i.e.,  the  rotational  axis  of  the  spheroid  is  parallel  to  the  direction  of  the  uniform  primary  field,  only  one  mode 
(0,  0, 1)  is  necessary  to  exactly  represent  the  field  as  shown  in  Fig.  2.4(a)  and  the  use  of  extra  terms  (p,  1, 1)  in 
Fig.  2.4(b)  or  even  higher  (not  shown  here)  has  no  contribution  to  this  decomposition  where  the  values  of  bpmn 
are  all  identically  zero.  When  the  spheroid  is  oriented  in  (ct,  9)  =  (0°,  45°),  the  two  modes  (0,  0, 1)  and  (0, 1, 1) 
are  needed  to  exactly  synthesize  the  field  //(,r .  In  this  case,  the  primary  field  equivalently  is  decomposed  into 
the  axial  and  transverse  excitations  in  the  body  coordinate-system.  The  field  behaves  like  a  two-dipole  model. 
For  the  spheroid  oriented  horizontally  along  x-axis,  i.e,  (a,  6)  =  (0°,90°),  this  is  a  configuration  of  transverse 
excitation.  Correspondingly,  we  see  that  in  Figs.  2.4(c)  and  (d)  the  contribution  is  exactly  zero  from  the  mode 
(0,0, 1)  and  only  the  mode  (0, 1, 1)  plays  a  role  in  the  decomposition.  The  primary  field  modes  (0,0, 1)  and 
(0, 1, 1)  correspond  to  axial  and  x-transverse  excitations,  respectively.  Similarly,  the  mode  (1, 1, 1)  corresponds 
to  y-transverse  excitation. 

With  a  uniform  field,  we  examine  its  decomposition  on  the  spheroid  oriented  in  different  angles.  The  exam¬ 
ples  are  simple  but  provide  some  physical  insight  into  the  mode  expansion  and  serve  to  check  the  stability  and 
correctness  of  the  algorithm  to  determine  the  bpmn. 
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(a)  Vertical:  0  =  0°  (b)  Vertical:  0  =  0° 


(c)  8  =  45° 


(d)  8  =  45° 


(e)  Horizontal  8  =  90°  (f)  Horizontal  8  =  90° 


Figure  2.4:  Spheroidal  expansion  of  a  uniform  field  H  =  z  on  the  spheroid  with  semi-major  and  minor  axes 
a  =  11  cm  and  6  =  5  cm.  The  orientation  of  the  spheroid  is  9  =  0°  for  (a)  and  (b);  9  =  45°  (c)  and  (d);  9  =  90° 
(e)  and  (f). 
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(a)  Vertical:  0  =  0° 


(b)  9  =  45° 


(c)  Horizontal  9  =  90° 


Figure  2.5:  Spheroidal  expansion  of  a  non-uniform  field  from  GEM-3  transmitter  on  the  spheroid  with  semi-major 
and  minor  axes  a  =  11  cm  and  b  =  5  cm.  The  vertical  distance  between  the  transmitter  and  the  center  of  the 
spheroid  is  29.2  cm.  The  orientation  of  the  spheroid  is  (a)  9  =  0°,  (b)  9  =  45°,  and  (c)  9  =  90°. 


2.4.2  Non-uniform  Primary  Field 

We  now  consider  the  expansion  of  a  spatially  non-uniform  field  produced  by  a  simulated  GEM-3  instrument.  The 
GEM-3  consists  of  two  concentric  circular  transmitter  coils  connected  in  series  but  with  opposite  current-flow 
direction.  By  properly  choosing  the  radius  of  the  two  coils  and  the  number  of  turns  of  the  coils  the  primary 
magnetic  flux  vanishes  in  the  receiver  loop.  Therefore,  the  small  receiving  coil  placed  within  this  magnetic  cavity 
senses  only  the  secondary  signal  returned  from  the  earth  and  nearby  metallic  objects. 

In  the  first  case,  we  consider  an  example  where  the  vertical  distance  between  the  transmitter  loops  and  the 
center  of  the  spheroid  is  29.2  cm.  Their  centers  are  aligned  in  the  same  vertical  line.  Fig.  2.5  presents  the  mode 
decomposition  for  the  spheroid  oriented  in  the  three  different  polar  angles.  More  terms  for  the  non-uniform  field 
are  necessary  to  accurately  express  the  field  than  the  case  for  the  uniform  field.  For  the  vertically  oriented  spheroid 
(i 9  =  0°),  the  accurate  expansion  requires  ( M,N )  =  (0,3).  For  the  spheroid  oriented  in  other  two  directions 
(0  =  45°), (9  =  90°),  we  need  the  maximum  number  (M,  N)  =  (1, 6)  and  (M,  N)  =  (4,  8)  to  almost  identically 
represent  the  field,  respectively.  For  the  transverse  case,  the  relative  error  is  around  0.33%  when  (M,  N)  =  (3,  8). 

When  the  vertical  distance  between  the  loop  and  the  center  of  the  spheroid  is  reduced  to  17  cm  (thereby  being 
the  same  scale  as  the  outer  loop),  the  degree  of  non-uniformity  is  increased.  We  see  in  Fig.  2.6  that  more  modes 
in  this  near  field  case  are  required  to  properly  represent  the  non-uniformity  of  the  near  field  pattern  produced  by 
the  GEM-3  transmitter.  For  a  transverse  excitation,  the  synthesized  fields  with  ( M ,  N )  =  (5, 8)  are  identical  with 
the  modeled  transmit  field  and  the  relative  error  is  around  0.40%  when  (M,  N)  =  (3,  8).  The  results  show  that 
more  N  modes  are  necessary  to  represent  a  non-uniform  field,  especially  for  a  transverse-like  excitation.  In  the 
following  application,  we  truncated  the  series  in  the  above  equations  at  (M,  N)  =  (3, 8). 
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-1  -0.5  0  0.5  1 

r)  coordinate  along  generating  arc 


(a)  Vertical:  0  =  0° 


(b)  9  =  45° 


(c)  Horizontal  9  =  90° 

Figure  2.6:  Spheroidal  expansion  of  a  non-uniform  field  from  GEM-3  transmitter  on  the  spheroid  with  semi-major 
and  minor  axes  a  =  11  cm  and  b  =  5  cm.  The  vertical  distance  between  the  transmitter  and  the  center  of  the 
spheroid  is  17  cm.  The  orientation  of  the  spheroid  is  (a)  8  =  0°,  (b)  8  =  45°,  and  (c)  8  =  90°. 

2.5  An  Example  of  SEA  Modeling 

Once  the  surface  sources  are  available  for  an  interested  target,  the  SEA  modeling  of  the  EMI  response  involves 
determining  the  spheroidal  modal  expansion  coefficients  bpmn  for  the  excitation  due  to  a  target-sensor  configura¬ 
tion.  We  demonstrate  this  by  modeling  the  RSS  for  a  solid  steel  cylinder  with  length  of  L  =  30.48  cm  and  the 
diameter  of  d  =  7.5  cm.  The  aspect  ratio  is  L/d  ss  4.06. 

We  consider  data  from  a  GEM3  frequency-domain  sensor  manufactured  by  Geophex  Ltd.,  and  time -domain 
data  collected  using  the  Geonics  Ltd.  EM63.  The  GEM3  data  were  collected  at  the  Sky  Research  Test  Plot  and 
the  EM63  data  were  collected  on  the  ERDC  test  stand. 

2.5.1  F  requency-domain  Response 

In  the  frequency-domain  measurements,  the  cylinder  was  excited  transversely  and  axially.  For  both  transverse  and 
axial  excitations,  the  vertical  distance  between  the  sensor  and  the  center  of  the  cylinder  is  h  =  42.75  cm  and  h  = 
45.24  cm,  respectively.  Figs.  2.7(a)  and  (b)  show  the  comparison  between  the  measured  and  modeled  data  with 
frequency  range  of  90  Hz  to  45  kHz  for  both  excitations.  The  SEA  modeling  was  very  successful  in  modeling  the 
measured  EMI  responses.  The  features  of  these  responses  conform  to  the  previous  studies  (Shubitidze  et  ah,  2002; 
Shubitidze,  O’Neill,  Sun  and  Paulsen,  2004).  The  crossover  frequencies  remain  essentially  fixed  in  Figs.  2.7(a) 
and  (b).  But  the  quadrature  peak  in  the  transverse  excitation  shifts  to  higher  frequency,  as  compared  to  that  in 
the  axial  excitation.  In  addition,  by  observing  the  quadrature  peak  frequencies  ft  and  fa  in  transverse  and  axial 
excitation,  we  can  estimate  the  aspect  ratio  as  \J  ft/  fa  =  i/l470/90  ~  4.04,  being  quite  close  to  the  actual  value 
of  4.06. 
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Figure  2.7:  Frequency-domain  response  from  GEM-3,  h  is  the  distance  between  the  center  of  the  cylinder  (L  = 
30.48  cm,  d  =  7.5  cm)  and  the  center  of  the  GEM-3  receiver. 


2.5.2  Transient  EM  Responses 

We  used  a  three  step  approach  to  modeling  the  TEM  response:  (1)  calculate  the  response  in  the  frequency  domain, 
(2)  use  a  Fourier  transform  to  convert  the  frequency  domain  response  to  a  time  domain  response  and,  (3)  convolve 
the  impulse  response  with  the  transmitter  waveform. 

The  step-off  and  impulse  responses  in  the  time  domain  can  be  evaluated  by  applying  the  digital  filter  tech¬ 
nique  (Anderson,  1982)  to  the  inverse  cosine  and  sine  transform  (Newman  et  al.,  1986)  as  follows 

2  r°°  l 

A(t)  — - /  ImB  (u)  —  cos  ojtduj ,  (2.9) 

7r  Jo  w 

where  A(t)  represents  a  step-off  response  corresponding  to  a  magnetic  flux  B(ui)  in  frequency  domain  and 
2  f°° 

A'(t)  =  —  /  ImB(w)  smutduj,  (2.10) 

7T  JO 

where  IuiB(lo)  represents  an  imaginary  part  of  the  magnetic  flux  and  the  prime  indicates  a  derivative  with  respect 
to  time  t. 

By  using  the  convolution  theorem  (Arfken,  1995),  the  voltage  in  the  time-domain  is 

^  =-J*  A' {t  -  r)I\r)dr  -  A(0)/'(f)  -  A'(t)I( 0).  (2.11) 

This  equation  represents  the  response  of  a  system  to  a  general  input  function  / (t )  in  terms  of  the  response  to  a 
unit  step  function  or  delta  function.  The  integration  in  Eq.  (2.1 1)  is  numerically  evaluated. 

The  transmitter  waveform  of  the  EM63  consists  of  an  exponential  current  increase  followed  by  a  linear  ramp 
off.  The  cylinder  was  oriented  at  polar  angles  of  6  =  0°,  45°,  90°.  For  each  of  the  three  excitations,  the  vertical 
distances  between  the  sensor  and  the  center  of  the  cylinder  are  h  —  60.00  cm.  Because  we  obtained  the  RSS  in 
a  wide  frequency  range,  the  transient  responses  were  produced  by  an  inverse  Fourier  transform  of  the  associated 
frequency-domain  signal  using  the  above  formulas.  The  comparison  in  Figs.  2.8(a)  -  (c)  shows  that  the  measured 
and  modeled  TEM  responses  agree  well  for  these  measured  data. 

2.6  Development  of  a  Reduced  Source  Library 

Once  we  determined  that  the  SEA  would  be  successful  in  modeling  the  sensor  response,  we  established  an  RSS- 
library  for  9  candidate  UXO,  as  listed  in  Table  2. 1  and  Fig. 2. 9.  These  items  include  UXO  from  the  Aberdeen  Test 
Center  (ATC)  UXO  test  set  and  items  provided  to  us  by  the  Montana  Army  National  Guard  (MTANG).  The  targets 
represent  a  mix  of  different  sized  ordnance  with  different  material  compositions  (steel,  aluminum,  copper),  and 
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(a)  Transverse  excitation 


E  10 


-  Measured 

-  -  SEA 


1  o~3  io-2 

Time  (sec) 

(b)  Axial  excitation 


(c)  Cylinder  with  a  dip  of  45  degrees 


Figure  2.8:  Transient  response  for  the  same  solid  steel  cylinder  in  Fig.  2.7.  For  each  measurement  the  cylinder  is 
60  cm  from  the  measurement  loop. 


include  projectiles,  submunitions  and  mortars.  Within  this  target  set,  UXO  (i.e.  a  2.75  inch  rocket  (labeled  C7) 
and  an  8 1  mm  mortar  (C8))  consisting  of  steel  and  aluminum  sections  would  be  challenging  to  model  accurately 
in  the  near  field  with  a  single  dipole  model. 


Candidate  No. 

Reference 

Material 

Total  length  (mm) 

Cl 

M55 

20  mm  steel  with  thin  cooper  band 

75 

C2 

M42 

40  mm  hollow  steel 

62 

C3 

M385 

40  mm  aluminium 

75 

C4 

M86 

57  mm  projectile,  thin  cooper  band  near  base 

170 

C5 

M49A3 

60  mm  steel 

130 

C6 

Montana 

76  mm  steel.  Armor  Piercing  with  Tracer  (APT) 

220 

C7 

M230 

2.75”  rocket,  aluminium  nose,  steel  body 

400 

C8 

M374 

81  mm  steel  body,  aluminium  nose  and  tail 

460 

C9 

Montana 

90  mm  projectile.  Target  Practice  with  Tracer  (TPT) 

250 

Table  2.1:  RSS-library  targets 

Each  target  was  photographed  and  digitized,  such  that  each  ordnance  could  be  modeled  as  a  body  of  revolu¬ 
tion  using  the  MAS.  Each  section  of  the  target  was  represented  by  a  cylinder.  Target  C7  and  C8  were  represented 
with  4  and  5  cylindrical  sections,  respectively.  The  conductivity  and  magnetic  permeability  of  each  section  were 
determined  by  fitting  data  from  a  GEM-3  sensor  collected  at  the  Sky  Research  Test  plot  in  Ashland,  July  2005 
(Figure  2.10).  Each  target  was  place  horizontally  beneath  the  GEM-3  sensor,  and  the  frequency  domain  response 
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(a)  Cl:  20  mm  M55 


(b)  C2:  M42  Submunition 


(c)  C3:  40  mm  M3  85 


(d)  C4:  57mm  M86 


(e)  C5:  60  mm  M49A3 


(f)  C6:  76  mm  from  Montana  Army  National 
Guard 


(g)  C7:  2.75  inch  Rocket,  M230  (h)  C8:  81  mmM374  (i)  C9:  90  mm  from  Montana  Army  National 

Guard 


Figure  2.9:  Items  for  which  RSS  were  computed. 
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Figure  2.10:  Apparatus  for  acquiring  data  for  determining  the  RSS. 


was  measured  at  operating  frequencies  of  90,  450,  930,  2010,  5970,  10050,  15450,  23010,  32010,  and  44010 
Hz.  Although  data  were  collected  along  a  line  over  the  target,  only  the  frequency  spectra  measured  directly  over 
the  center  of  the  target  was  used  for  the  fitting.  Fitting  these  data  was  achieved  through  ’’trial  and  error”.  The 
off-center  data  were  used  as  a  check  to  determine  the  suitability  of  the  selected  permeability  and  susceptibility  for 
modeling.  Clearly,  this  procedure  needs  to  be  refined  to  ensure  greater  accuracy  and  repeatability. 

The  RSS  library  was  built  using  three  different  sizes  of  spheroids  with  a  =  11.0  cm,  b  =  5.0  cm  for  small 
objects,  a  =  22.0  cm,  b  =  7.5  cm  for  medium  ones,  and  a  =  30.0  cm,  b  =  7.5  cm  for  the  largest  items.  On 
the  surface  of  these  spheroids,  12  rings  of  charge  were  used.  Sun  et  al.  (2005)  used  6  rings  in  their  analysis. 
Figures  2.11  and  2.12  contain  the  RSS  models  as  a  function  of  length  along  the  spheroidal  coordinate  //,  for 
targets  C3  (40  mm  aluminum  submution)  and  C5  (60  mm  steel  mortar  body)  respectively. 

Figure  2.13  compares  the  response  of  an  40  mm  (C3)  predicted  using  the  RSS  as  depicted  in  Figure  2.11 
assuming  a  GEM3  sensor,  with  data  collected  on  the  ERDC  test  stand.  The  target  depth  was  30  cm  from  the  plane 
on  which  the  GEM3  head  moves.  This  target  represents  the  most  straightforward  object  from  our  library  to  model. 
It  does  not  contain  a  number  of  different  sections  and  is  non-magnetic.  The  plan  view  images  of  Figure  2.13(a) 
and  (b)  show  the  real  and  imaginary  parts  for  3  of  the  10  frequencies.  The  line  plot  at  the  bottom  of  each  panel 
plots  the  frequency  response  measured  and  modeled  directly  above  the  center  of  the  target  (i.e.  (X,  Y)  =  (0, 0). 
This  target  response  was  very  accurately  modeled  using  the  RSS  from  our  library. 

Figure  2.14  compares  the  response  of  an  81  mm  (C8)  modeled  with  the  RSS  and  GEM3  sensor  data.  The 
target  depth  was  50  cm  from  the  plane  on  which  the  GEM3  moves.  This  target  represents  a  much  more  challenging 
item  to  model.  The  ATC  81  mm  contained  several  different  sections;  both  ferrous  and  non-ferrous. 

The  crossover  of  the  real  and  imaginary  parts  is  slightly  different  between  measured  response  and  modeled 
response.  However,  as  is  seen  in  the  plan  data,  at  the  higher  and  lower  frequencies  the  plan  view  images,  and 
therefore  the  data,  match  closely.  The  slight  inaccuracy  of  the  crossover  suggests  more  care  should  be  made 
when  calculating  the  RSS  for  the  larger,  more  complex  targets.  Regardless,  the  comparison  between  predicted 
and  measured  data  is  quite  good,  and  we  proceed  with  testing  discrimination  using  our  RSS  library. 
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RSS  amplitude  RSS  amplitude  RSS  amplitude  RSS  amplitude 


(d)  /  =  10481.13  Hz:  m  =  0,  n  =  1 


(e)  /  =  10481.13  Hz:  m  =  1,  n  =  1 


(f)  /  =  10481.13  Hz:  m  =  3,  n  =  4 


Figure  2.1 1:  The  distribution  of  RSS  for  C3  along  a  spheroidal  coordinate  rj. 


(a)  /  =  126.48  Hz:  m  =  0,  n  =  1 


(b)  /  =  126.48  Hz:  m  =  1.  n.  =  1 


(c)  /  =  126.48  Hz:  m  =  3,  n  =  5 


Figure  2.12:  The  distribution  of  RSS  for  C5  along  a  spheroidal  coordinate  rj. 
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(a)  RSS  generated  data 


(b)  Data  measurd  by  the  GEM3  sensor 
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Figure  2.13:  Modeled  and  measured  GEM3  response  for  a  horizontal  ATC  40  mm  (C3). 


(b)  Data  measurd  by  the  GEM3  sensor 


(a)  RSS  generated  data 
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Figure  2.14:  Modeled  and  measured  GEM3  response  for  a  horizontally  ATC  81mm  (C8). 


2.7  UXO  Identification  Tests 


The  above  examples  illustrate  the  applicability  of  the  MAS/TSA-based  SEA  for  the  fast  and  accurate  EMI  mod¬ 
eling.  This  modeling  technique  is  readily  inserted  into  ’’fingerprinting”  or  template  matching  type  discrimina¬ 
tion/classification  algorithms.  The  objective  of  our  template  matching  analysis  is  to  determine,  from  a  list  of  M 
targets,  the  target  that  is  most  likely  to  have  generated  the  observed  data  dobs.  Each  target  in  our  list  is  charac¬ 
terized  by  its  RSS,  which  is  represented  by  the  vector  p,  .  For  each  RSS  in  our  library,  we  determine  the  location 
r.i  and  orientation  angles  0,  and  (Ip,  at  which  we  can  best  fit  the  observed  data  dobs  by  obtaining  the  maximum 
likelihood  solution.  The  data  predicted  by  this  recovered  model,  dfred  =  F[r,,  0, .  <pt.  p,]  =  is  referred 

to  as  the  template  for  target  i.  The  target  template  dfred  most  similar  to  the  observed  data  dobs  is  selected  as  the 
most  likely  target. 

There  are  a  number  of  measures  that  compare  the  target  templates  with  the  observed  data.  Intuitively,  these 
include  measures  of  maximum  correlation  or  minimum  error.  There  are  also  several  ways  with  which  to  define 
the  minimum  error.  Riggs  et  al.  (2001)  outline  the  derivation  of  the  minimum  least  squares  from  a  generalized 
likelihood  ratio  test  (GLRT)  with  Gaussian  data  statistics.  The  likelihood  ratio  test  for  two  targets 

p(dobs|pi)  t°.rgetl  p(p2)  (C01  -  Cn)  =  ?  ^ 

V  (dobs|p2)  target2  V  (Pi)  (C12  —  C22) 

where  C,:/  is  the  cost  of  classifying  the  target  as  p,  when  the  target  is  p(,  and  p  (p,  )  is  the  prior  probability 
for  the  ith  class.  The  GLRT  is  obtained  by  substituting  the  maximum  likelihood  estimate  into  equation  2.12. 
If  we  consider  two  targets  with  equal  prior  probability  of  producing  the  anomaly,  and  assuming  that  an  incorrect 
classification  produces  the  same  cost,  then  77  =  1.  By  taking  the  logarithm  of  the  resulting  expression,  our  decision 
criterion  is  to  simply  select  the  target  that  has  the  smallest  least  squares  error: 

||k71/2(dobs-F[r1,01,  ^Pi])  ||2  ^  ||k71/2(dobs-  F[r2,02,<h,p2)])  ||2  (2.13) 

target2 

where  r,  is  the  position,  and  (/,.  and  (]>t  are  the  orientation  angles  that  produces  the  best  fit  to  the  observed  data  for 
the  model  p,  .  For  multiple  candidate  targets  we  simply  choose  the  target  with  smallest  least  squares  misfit. 

Application  of  the  above  algorithm  to  survey  data  requires  establishing  two  thresholds.  First,  a  minimum 
level  of  data  quality  must  be  established,  since  our  confidence  in  identifying  the  correct  target  decreases  as  the  data 
quality  degrades.  For  data  sets  acquired  with  the  same  survey  parameters  (such  as  station  density)  for  the  entire 
data  set,  the  critical  measure  of  data  quality  is  the  signal  to  noise  ratio  (SNR).  A  second  threshold  is  a  maximum 
misfit  at  which  an  anomaly  can  be  labeled  as  a  target  within  our  library.  Since  the  ability  to  distinguish  differences 
between  the  observed  data  and  the  template  data  will  depend  on  the  quality  of  data,  the  minimum  correlation 
threshold  will  also  be  dependent  on  the  Signal-to-Noise  Ratio  (SNR)  ratio  of  the  target.  These  thresholds  can  be 
established  with  training  data  or,  if  the  survey  noise  can  be  accurately  modeled,  through  simulations. 


2.7.1  A  simple  discrimination  example  using  a  single  sounding 

For  a  first  test,  we  consider  a  very  simplified  example.  For  this  test  we  use  the  9  candidate  patterns  of  the  scattered 
fields  measured  at  different  levels  ranging  from  15  cm  to  30  cm  above  the  objects,  being  positioned  at  (0,  0,  0) 
cm  and  oriented  at  (a,  9)  =  (0°,  90°).  This  example  does  not  represent  a  realistic  discrimination  test,  but  rather 
serves  as  a  simple  test  of  the  variations  in  soundings  that  are  generated  by  the  different  members  of  the  library 
and  illustrates  the  basic  idea  of  the  pattern  matching.  To  be  consistent  with  Sun  et  al.  (2005),  we  characterize  the 
difference  between  the  trial  and  candidate  patterns  with  the  following  misfit  function 


F{Ci,Cj) 


E/l^obs(/,C'i)-^trial(/,a)| 

£/  \Hohs(f,Ci)\ 


(2.14) 


where  fT0bs(/:  Ci)  are  the  measured  candidate  pattern  of  the  scattered  fields  for  unknown  candidate  C,  to  be 
identified,  and  //trial  (/,  C:l  )  are  the  trial  pattern  of  the  secondary  fields  for  UXO  Cj  on  the  basis  of  the  RSS- 
library. 

Given  a  candidate  pattern,  for  instance,  we  can  quickly  scan  the  trial  patterns  with  it  and  produce  the  misfit 
values  as  shown  in  Fig.  2.15.  It  can  be  observed  that  the  minimum  misfit  value  occurs  for  candidate  C4.  In  fact 
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Canadiate  No. 

Figure  2.15:  Misfit  values  for  a  trial  identification  test. 


the  value  of  F(Ci ,  C4)  is  only  2.2%;  this  is  far  smaller  than  other  misfit  values  of  F(Ci ,  Cj),j  ^  4.  Assuming 
an  error  threshold  of  10%,  we  are  able  to  identify  this  unknown  object  as  UXO  C4  corresponding  to  the  given 
candidate  pattern.  This  identification  result  is  confirmed  by  visually  inspecting  the  patterns  of  the  secondary  fields 
between  the  calculated  and  the  candidate  in  Fig.  2.16. 

In  another  trial  test,  we  sweep  the  other  candidate  pattern  using  the  same  RSS-library.  The  misfit  values  are 
presented  in  Fig.  2.17.  Again,  one  can  readily  locate  a  minimum  value  of  F(Ci,  C's )  -  being  3.4%  below  the  error 
threshold  value.  The  other  values  of  F{Ci ,  Cj),j  7^  8  are  much  larger.  Based  on  this,  we  can  make  a  decision  to 
identify  this  unknown  object  as  UXO  C8.  Fig.  2.18  illustrates  that  UXO  C8  is  the  only  item  with  a  good  match 
between  the  trial  and  the  measured  data. 
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Figure  2.16:  Pattern  matching  tests,  (a)-(i):  C1-C9.  UXO  C4  identified 
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Figure  2.17:  Misfit  values  for  another  trial  identification  test. 
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Figure  2.18:  Pattern  matching  tests.  UXO  C8  identified 
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2.7.2  A  discrimination  example  using  GEM3  data  collected  on  the  ERDC  Test  Stand 

The  example  presented  in  the  previous  section  is  not  realistic,  since  location  and  orientation  are  not  generally 
known  prior  to  processing  data.  Although  a  priori  information  on  position  and  orientation  can  sometimes  be 
obtained  from  prior  surveys,  often  the  information  is  inaccurate.  The  initial  plan  was  to  implement  the  same 
library  based  identification  techniques  as  in  Pasion  et  al.  (2006).  In  that  approach,  a  library  of  polarization  tensors 
was  created  from  time  domain  electromagnetic  data  collected  on  a  test  stand.  For  each  target  in  our  library  a 
non-linear  inverse  problem  was  solved  for  the  position  and  orientation  that  minimized  the  least-squares  difference 
between  the  observed  data  anomaly  and  the  data  predicted  from  each  target.  Our  code  for  inverting  for  location 
and  orientation  is  being  developed  and  is  not  mature  enough  to  implement  the  non-linear  inversion  approach. 
On  the  other  hand,  our  initially  optimized  RSS  forward  modeling  code  is  able  to  rapidly  yield  the  response  with 
around  2.86  seconds  per  sounding  of  13  frequencies.  Therefore,  we  developed  a  template  matching  technique  to 
determine  depth  and  orientation  by  searching  a  library  of  data  pre-modeled  at  several  depths  and  dip  angles  and 
used  image  registration  to  find  the  location  and  azimuth  angle  of  the  target.  The  main  objective  of  implementing 
this  style  of  template  matching  is  to  determine  if  it  is  possible  to  identify  targets  using  an  RSS  library  without  a 
priori  information.  We  emphasize  that  we  do  not  claim  the  algorithm  presented  here  is  the  best  choice  for  target 
identification,  and  we  plan  to  implement  the  non-linear  inverse  problem  approach  in  Pasion  et  al.  (2006)  at  a  later 
date. 

The  first  step  of  the  algorithm  is  to  generate  a  library  of  UXO  responses.  The  generation  of  a  library  of 
UXO  responses  meant  that  all  the  forward  modeling  using  the  RSS  only  needed  to  be  performed  once,  thereby 
increasing  the  speed  of  the  analysis.  Target  responses  for  items  Cl  to  C9  were  calculated  for  target  distances  from 
the  GEM3  sensor  head  varying,  at  10  cm  intervals,  from  20  cm  to  80  cm.  At  each  depth  the  target  was  measured 
at  dip  angles  from  0  degrees  (horizontal)  to  90  degrees  (vertical),  at  15  degree  intervals.  Data  were  modeled  on  1 
m  square  area  and  on  a  uniform  grid  with  10  cm  spacing. 

UXO  identification  is  achieved  by  determining  the  data  template  from  our  library  that  best  matches  the 
sensor  data  by  cycling  through  each  of  the  data  templates.  However,  the  target  location  and  azimuthal  orientation 
is  unknown.  Determining  the  target  location  and  azimuthal  orientation  is  equivalent  to  determining  the  translation 
and  rotation  of  the  data  templates.  This  operation  represents  a  simple  problem  in  image  registration,  since  we  do 
not  have  to  consider  scaling  the  template.  Figure  2.19  shows  an  example  of  determining  the  correct  translation, 
and  therefore  an  estimate  of  location.  The  upper  left  panel  shows  sensor  data  collected  on  a  test  stand.  The 
upper  right  panel  shows  the  1  m  square  area  of  the  template.  Both  the  data  and  the  template  are  gridded  to  2.5 
cm  pixels.  The  lower  left  panel  plots  the  correlation  coefficient  as  a  function  of  position.  Not  surprisingly,  the 
optimal  location  is  directly  over  the  center.  In  the  examples  that  follow,  we  do  not  implement  an  image  registration 
technique  and  we  assume  that  we  know  the  azimuthal  orientation  and  location.  We  note  that  while  implementing 
an  image  registration  algorithm  is  straightforward,  we  would  like  to  focus  on  our  ability  to  correctly  determine 
depth  and  dip. 

To  demonstrate  the  above  procedure,  we  used  GEM3  data  collected  on  the  USACE-ERDC  Test  Stand.  The 
data  were  collected  from  January  31,  2006  to  February  10,  2006.  For  a  first  example  we  compare  data  from  a 
vertical  40  mm  projectile  located  30  cm  from  the  GEM3  sensor  head.  Figure  2.20  compares  images  of  measured 
sensor  data  with  the  best  fit  template  for  3  of  the  10  GEM3  frequencies.  Figure  2.21(a)  compares  the  sounding 
directly  over  the  center  of  the  target.  Figure  2.21(b)  compares  the  misfit  values  for  the  different  items  in  the  RSS 
library.  It  is  clear  that  C3  is  the  most  likely  target.  The  ease  with  which  the  algorithm  picked  C3  is,  in  part,  due  to 
C3  being  the  only  non-ferrous  (aluminum)  item  in  the  library. 

Figures  2.22  and  2.23  plots  the  results  when  the  data  from  an  81  mm  (M374)  is  fit.  Although  the  correct 
target  (C8)  has  the  minimum  misfit,  several  other  targets  have  relatively  similar  misfit  values,  indicating  that 
several  targets  appear  quite  similar  when  viewed  by  a  GEM3  sensor  50  cm  from  the  target. 
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Figure  2.19:  Examples  of  finding  the  optimal  location  by  translating  a  template  generated  by  the  RSS. 
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Figure  2.20:  Modeled  and  measured  GEM3  response  for  a  vertically  oriented  ATC  40  mm  (C3) 
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Figure  2.21:  Template  example  for  C3. 
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Figure  2.22:  Modeled  and  measured  GEM3  response  for  a  horizontally  oriented  ATC  81  mm  (C8). 
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Figure  2.23:  Template  example  for  C8. 


2.8  Inversion  for  Optimal  Orientation  and  Location  of  an  Object 

The  above  identification  tests  were  conducted  assuming  an  object  location  and  orientation  as  a  priori  information. 
When  such  information  is  unavailable  such  as  in  a  usual  discrimination  scenario,  we  need  to  estimate  them  from 
sensor  measurements.  This  section  reports  such  a  preliminary  test  to  invert  the  optimal  location  and  orientation  of 
a  known  object  by  using  a  gradient-type  method. 

Denote  the  parameter  vector  m  =  (rc,  a,  9),  where  rc  denotes  the  position  vector  of  an  object  center,  a 
and  9  are  the  azimuthal  and  polar  angles  of  the  object  relative  to  the  user  coordinate.  Based  on  the  Eq.  (2.4) 
of  the  secondary  magnetic  field,  we  can  derive  the  analytical  expression  for  the  sensitivities  of  Hsc  to  the  5 
parameters.  For  this  nonlinear  inverse  problem,  a  optimization  approach  is  implemented  iteratively  to  correct  the 
model  parameters  via  minimizing  the  objective  functional 

F( TC,  a ,  9)  =  ||HPcre(rc,  a,  9)  -  H°cbs||2,  (2.15) 

where  TF’”’  are  the  predicted  fields  at  a  trial  model,  H°bs  are  the  observed  ones. 

The  cylinder  in  section  4  is  chosen  as  an  object,  being  positioned  at  (0, 0,  —42.5)  cm  and  oriented  at  9  =  7t/4 
and  a  =  7t/4.  For  this  inversion,  the  7x7  simulated  GEM-3  measurement  points  are  distributed  on  a  plane  grid 
with  the  interval  of  20  cm  and  a  single  frequency  of  41010  Hz  is  used  in  the  inversion.  In  the  inversion  test,  we  put 
the  lower  and  upper  bound  vectors  m;  =  [—30  cm  —  30  cm  —  70  cm  0  0]  and  mu  =  [30  cm  30  cm  —  30  cm  27t  zt]  . 
Fig.  2.24  (a)  -  (f)  illustrates  the  convergence  of  the  inverted  results  versus  the  iteration  number.  It  is  observed 
that  the  after  10  iterations,  the  model  parameters  are  recovered  very  well.  The  values  of  data  misfit  monotonically 
decrease  with  the  iteration  and  have  no  significant  variations  after  10  iterations,  as  shown  in  Fig.  2.24  (f). 

In  another  test,  the  cylinder  takes  the  same  position  as  before  but  oriented  at  9  =  0°  and  a  =  0°.  Fig. 
2.25  shows  inverted  results  that  are  fairly  good.  The  convergence  behavior  in  Fig.  2.25  (f)  is  similar  to  that  of 
the  previous  test.  Note  that  the  recovered  value  of  9  finally  approaches  to  tt  not  0  because  both  values  have  no 
difference  for  a  homogeneous  cylinder.  In  this  case  the  value  of  the  azimuthal  angle  a  could  be  arbitrary  as  the 
object  is  oriented  vertically. 

2.9  Conclusions 

We  have  investigated  the  MAS/TSA-based  implementation  of  the  standardized  (spheroidal)  excitation  approach 
for  a  full  3D  electromagnetic  induction  sensing  of  metallic  objects.  The  approach,  by  decomposing  excitation 
fields  into  fundamental  spheroidal  modes,  synthesizes  the  secondary  fields  through  a  linear  superposition  of  each 
modal  response  to  a  target.  We  quantitatively  studied  the  fundamental  issue  of  determining  the  effective  spheroids 
in  decomposition.  It  was  found  that  for  a  non-uniform  primary  field  excited  by  a  GEM3-like  loop,  the  infinite 
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Figure  2.24:  Inversion  test  I.  Recovered  model  parameters  versus  iteration  number. 


series  can  be  truncated  at  the  maximum  mode  number  (M,  N)  =  (3, 8),  which  is  used  to  accurately  represent  the 
fields  for  various  orientations  of  small/medium  objects  in  a  near  field  zone. 

The  key  feature  in  this  approach  is  that  the  extrinsic  and  intrinsic  attributes  characterizing  the  pattern  of 
secondary  fields  are  automatically  separated  out  in  the  spheroidal  mode  expansion  coefficients  and  the  modal 
response  amplitudes  or  the  RSS,  respectively.  Based  on  the  RSS  for  interested  objects,  the  related  full  3D  EMI 
modeling  can  be  carried  out  rapidly  by  merely  determining  primary  field  expansion  coefficients  for  any  excitation- 
object  geometry.  These  were  illustrated  by  the  tests  against  the  experimental  data  in  both  frequency-domain 
(GEM-3)  and  time-domain  (EM63)  systems.  We  extended  the  existing  SEA  to  simulate  a  transient  EM  response 
by  developing  a  convolution  algorithm  that  applies  to  an  arbitrary  current  waveform  used  in  a  TEM  system. 
These  experiment-based  results  are  important  not  only  to  confirm  the  soundness  of  the  approach  but  also  indicate 
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Figure  2.25:  Inversion  test  II.  Recovered  model  parameters  versus  iteration  number. 


its  promise  in  real-world  applications. 

Since  the  MAS/TSA-based  SEA  is  a  physically  complete  numerical  technique,  as  stated  previously,  it  can 
be  used  to  accurately  build  up  a  RSS-library  in  an  ultra  wide  band  -  ranging  from  resistive  to  inductive  limit.  With 
such  an  ultra  wide  band  RSS-library,  it  is  convenient  to  compute  EMI  response  in  either  the  time  or  frequency 
domain. 

As  another  direct  application  of  a  RSS-library,  we  have  conducted  identification  test  of  UXO  based  on  the 
test  stand  data.  In  this  test,  we  successfully  showed  that  identifying  a  candidate  UXO  can  be  quickly  done  through 
a  straightforward  pattern  matching  procedure  by  inspecting  a  best  misfit  values  below  an  error  threshold.  In  our 
follow-up  research  we  will  proceed  to  do  a  blind  test  of  the  approach  combined  with  the  robust  inversion  algorithm 
determining  an  optimal  object  location  and  orientation. 
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To  conclude,  the  MAS/TSA-based  SEA  is  a  promising  full  3D  EMI  modeling  technique  that  provides  a 
platform  on  which  to  develop  discrimination  and  identification  algorithms. 
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Chapter  3 

Regularization  of  the  Surface  Charge 
Model  for  the  Inversion  of 
Electromagnetic  Data 

3.1  Abstract 


Detection  of  buried  unexploded  ordnance  with  electromagnetic  sensors  requires  robust  predictive  models  to  cor¬ 
rectly  interpret  data  recorded  at  the  surface.  Shubitidze  et  al.  (2005a)  have  recently  suggested  a  physics-based 
representation  of  the  scattered  electromagnetic  field  using  surface  magnetic  charges.  One  such  representation 
is  the  Normalized  Surface  Magnetic  Charge  (NSMC)  model,  which  idealizes  a  compact  metallic  target  with 
magnetic  charges  distributed  on  a  spheroidal  surface  that  encloses  the  target.  The  Total  Normalized  Magnetic 
Charge  (TNMC),  as  a  function  of  frequency  or  as  a  function  of  time,  can  be  used  to  identify  the  object.  At  a 
theoretical  level,  however,  the  inference  of  the  TNMC,  is  an  ill-determined  problem.  Some  form  of  regularization 
of  the  problem  is  required.  In  order  to  help  develop  the  NSMC  model  into  a  robust  practical  tool  for  UXO  discrim¬ 
ination,  we  conducted  a  theoretical  study  on  the  behavior  of  these  magnetic  surface  charges  and  formulated  the 
calculation  of  the  total  normalized  charge  as  an  inverse  problem.  Extensive  use  was  made  of  Geophex  GEM3  and 
Geonics  EM63  test  stand  data  collected  at  the  USACE-ERDC  test  stand  in  Vicksburg,  MS.  Results  are  presented 
for  both  test  stand  data  and  data  collected  on  the  Sky  Research  UXO  Test  Site  in  Ashland,  Oregon. 


3.2  Introduction 


In  this  chapter  we  consider  the  Surface  Magnetic  Charge  (SMC)  model  of  Shubitidze  et  al.  (2005a).  The  SMC 
model  approximates  a  target’s  secondary  electromagnetic  field  at  each  frequency  or  time  with  the  field  from  a 
magnetic  charge  distribution  on  a  fictitious  surface  that  encloses  the  target  (Figure  3.1).  Shubitidze  et  al.  (2005a) 
claims  that  normalizing  the  charge  distribution  by  the  normal  component  of  the  primary  field  produces  a  surface 
distribution  that  is  a  property  of  the  target,  and  not  the  excitation  of  the  target.  Furthermore,  the  surface  integral 
of  the  normalized  charge  distribution  was  suggested  as  a  discriminant  for  UXO. 

Unlike  the  SEA  model  of  the  previous  chapter,  the  nature  of  SMC  and  its  computational  speed  make  it  a 
candidate  for  data  inversion  for  the  charge  distribution.  In  this  work,  we  study  some  of  the  properties  of  the 
magnetic  charge  model  and  some  of  the  properties  of  the  SMC  inverse  problem.  The  objective  of  this  work  was 
to  determine  if  a  stable  discriminant  derived  from  the  surface  magnetic  charge  model  could  be  determined. 


36 


(a)  The  dipole  model,  showing  axial  and  transverse  polarizations  (b)  The  Surface  Magnetic  Charge  (SMC)  model 


Figure  3.1:  Representation  of  the  response  of  a  target  using  (a)  a  dipole  and  (b)  a  surface  magnetic  charge. 

3.3  Model  description  and  implementation 

3.3.1  The  Surface  Magnetic  Charge  Model 


The  following  description  of  the  SMC  model  follows  Shubitidze  et  al.  (2005a).  They  assume  a  highly  conducting, 
magnetically  permeable,  heterogeneous  metallic  target  buried  in  soil  where  conductivity  is  low  enough  to  be 
considered  as  free  space.  In  a  quasi-magneto-static  regime,  displacement  currents  are  negligible,  conduction 
currents  are  weak  outside  the  target,  the  magnetic  field  is  irrotational  and  can  be  written  as  the  gradient  of  a  scalar 
potential  %l>\ 

Hsc(r,£)  =  -Vvrc(r,£)  (3.1) 

where  the  variable  £  can  represent  either  time  t  or  frequency  to.  Assuming  Gauss’  Law  for  magnetic  field  yields 
V  •  MoH  =  V2T-sc  =  pm  (3.2) 


where  pm  represents  a  fictitious  volumetric  magnetic  charge  density  that  produces  a  scattered  field  outside  a 
metallic  body.  The  scalar  potential  satisfies  an  equation  of  the  type 


*sc(r,0 


1 

47 xpo 


Pm  (r',0 


r  —  r 


rdV' 


(3.3) 


where  r  is  the  observation  point,  r'  the  source  point,  V  the  volume  of  the  scatterer,  po  the  magnetic  permeability 
of  free  space.  Assuming  surface  charges  am  only,  the  magnetic  field  is  given  by 


Hsc  (r,  £) 


1 

Airpo 


dS’ 


(3.4) 


The  charge  distribution  am  is  a  function  of  the  field  at  the  surface  S  which  is,  in  turn,  a  function  of  the  target 
properties  and  primary  field. 


The  surface  charge  is  assumed  to  be  proportional  to  the  normal  component  of  the  incident  primary  magnetic 

field: 

Cm  (r',  £)  =  qm  O',  0  [Hpr  O')  •  n  (r')]  (3.5) 

where  qm  is  the  NSMC  density  that  we  try  to  calculate.  Considering  that  UXO  are  bodies  of  revolution  with 
symmetry  along  their  main  axis,  Shubitidze  et  al.  (2005b)  assume  that  (1)  the  electromagnetic  response  of  a  UXO 
can  be  equivalently  represented  by  surface  charges  located  on  a  prolate  spheroid  enclosing  the  ordnance  with  same 
main  axis,  and  (2)  that  qm  on  this  spheroid  is  uniform  on  rings  perpendicular  to  the  main  axis.  The  TNMC  of  an 
object  at  a  given  time  or  frequency  is  defined  as 

Q{0=  f  qm(r',0  dS'  (3.6) 

Js 

Shubitidze  et  al.  (2005a)  suggested  that  the  TNMC  is  unique  for  each  target  so  it  can  be  used  to  discriminate 
between  different  targets. 
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3.3.2  Numerical  Implementation  of  the  SMC  Model 


For  this  report,  we  only  consider  fields  measured  by  induction  loop  receivers.  For  a  frequency  domain  system,  the 
voltage  induced  in  a  receiver  loop  is 

V  M  =  iwnR<f>  (3.7) 


where  nR  is  the  number  of  turns  in  the  receiver,  $  is  the  flux  through  the  receiver.  The  corresponding  expression 
in  the  time  domain  is 

<9<f> 

V{t)=nR—  (3.8) 

at 

The  flux  through  a  receiver  loop  in  air  is 


$  (r,  0  =  Ho  f  (Hsc  (0  •  np)  dSR  (3.9) 

JsR 

where  nR  is  the  normal  of  the  loop  surface,  and  we  again  define  £  as  being  either  t  or  u.  Assuming  a  horizontal 
loop: 

(310) 

=  ±f  I  0’  r')  0  -  z')  9  (r':  0  [Hp  •  n  (r')]  dS' 

where 

I(  r,r')=  [  —i—^dSR  (3.11) 

Js*  I r  —  r'  | 3 

depends  on  the  shape  of  the  loop.  We  solve  this  integral  analytically  for  a  horizontal,  rectangular  loop.  For  a 
horizontal  circular  loop,  we  solve  the  integral  in  polar  coordinates,  with  an  analytic  integration  over  the  radial 
coordinate,  and  a  numerical  integration  over  the  azimuthal  coordinate. 


Equation  3.10  shows  that,  once  the  charge  surface  is  determined,  the  flux  <f>,  and  therefore  the  induced 
voltage  V  is  a  linear  function  of  the  charge  distribution  q.  To  numerically  solve  the  SMC  integral,  we  discretize 
the  surface  S  into  N  sub-surfaces,  with  area  AS),  and  assume  that  the  amplitude  of  the  surface  magnetic  charge 
qi  is  constant  in  A s».  Substitution  in  Equation  3.10  yields 

1  N 

$(r’£)  =  ^  (z  -  zk)<l(rk,  0  [Hp  (rfc)  ■  n(rfc)]  A Sk  (3.12) 

k 

The  forward  model  data  collected  at  N  stations  at  a  frequency  uji  can  be  written  as 

Z  m  =  d  (3.13) 

where 

[z]mn  (r"»rn)  (zrn  -  zn)  [HP  (rn)  •  n  (rn)]  A Sn  (3.14) 

m  =  [q(r1,wi)  ,  q(r2,Ui)  ,  ,  q(rN,Ui)]T  (3.15) 

and  an  element  of  the  data  vector  d  is  the  integrated  secondary  field  flux  over  the  receiver  loop. 


For  time  domain  data. 


dq  (ri,  fj)  dq(ri,ti)  dq(TUti)' 

m "  [  at  ’  at  ’  '  ■  ‘  ’  at 

d  =  [V  (ri ,U)  ,  V  (r2, U)  ,  . . .  ,  V  (r N,ti)]T  . 


(3.16) 

(3.17) 
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The  modeling  matrix  Z  is  a  function  of  the  shape  of  the  charge  surface.  For  a  spheroid  that  contains  the 
axi-symmetric  target,  Z  is  then  a  function  of  the  location  and  orientation.  In  the  case  of  a  sphere  Z  is  a  function 
of  position  only. 


3.3.3  Discretization  of  the  charge  surface 


In  order  to  solve  for  a  representative  charge  distribution,  a  surface  must  be  defined  and  discretized.  In  this  study 
we  chose  spherical  and  spheroidal  surfaces  to  enclose  axi-symmetric  targets.  Figure  3.2  illustrates  a  discretization 
on  a  spheroid.  Patches  are  defined  by  uniformly  discretizing  in  the  azimuthal  coordinate  (longitude)  and  in  the 
local  axial  (z)  coordinate  (or  in  the  latitude  angle).  Rings  are  defined  by  assuming  that  the  surface  charge  density 
remains  an  invariant  of  the  azimuth. 


Figure  3.2:  An  example  of  discretizing  a  spheroid. 

Discretization  along  longitudinal  and  latitudinal  directions  samples  the  surface  of  a  spheroid  in  a  non- 
uniform  manner,  resulting  in  high  density  of  surface  patches  at  the  poles.  To  circumvent  that  potential  source 
of  error  we  tested  the  SMC  method  on  a  sphere  by  using  a  geodesic  discretization  derived  from  refinements  of  an 
icosahedron  (polyhedron  with  12  vertices,  20  triangular  surfaces  of  equal  area,  30  edges),  illustrated  in  Figure  3.3, 
where  red  lines  show  the  projection  of  the  original  icosahedron  onto  its  surrounding  sphere.  This  configuration  is 
almost  isotropic  and  every  patch  has  the  exact  same  dimensions,  thus  avoiding  any  area-induced  bias  in  the  charge 
density. 


3.4  Inverting  for  a  Surface  Charge  Distribution 


In  order  for  the  SMC  formulation  to  be  part  of  a  successful  discrimination  algorithm,  the  recovered  surface 
magnetic  charge  model  (or  some  property  of  the  SMC)  must  be  a  function  of  the  target  characteristics  only,  and 
independent  of  the  survey  parameters.  We  studied  the  sensitivity  of  the  SMC  solution  to  different  modeling 
parameters  such  as  depth,  orientation  and  noise  by  using  test  stand  data  and  synthetic  data  obtained  with  a  dipole 
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Refinement  T2 


(a)  80  faces  (b)  1280  faces 

Figure  3.3:  Discretization  of  a  sphere  derived  from  an  icosahedron 


model.  We  find  that  for  the  type  of  objects  of  interest,  and  for  the  type  of  data  that  are  collected,  the  problem 
of  solving  the  equation  Zm  =  d  is  often  ill-posed,  even  in  over-determined  cases  where  there  are  more  data 
than  parameters,  as  is  the  situation  with  the  SMC  method  of  rings.  In  the  more  general  under-determined  case, 
where  patches  cover  the  spheroid,  there  are  more  parameters  than  data  and  the  inverse  problem  is  severely  ill- 
posed.  Thus  there  are  infinitely  many  solutions  for  the  matrix  system  and  some  additional  information  needs  to 
be  incorporated  to  solve  the  problem.  The  instability  of  the  inverse  problem,  and  the  non-uniqueness  of  solution, 
is  exemplified  by  the  condition  number  of  Z,  cond  (Z),  the  ratio  of  the  largest  and  smallest  non-zero  singular 
values  of  Z.  The  condition  number  grows  rapidly  with  the  resolution  of  the  model,  for  instance  with  10  rings 
cond  (Z)  =  10s,  with  40  rings  cond  (Z)  =  101'  and  with  200  patches  cond  (Z)  =  1013. 

We  obtained  a  least  squares  solution  to  the  matrix  equation  by  using  the  pinv  command  in  Matlab,  which 
computes  the  Moore-Penrose  pseudo-inverse  of  a  matrix,  effectively  solving  ZTZm  =  ZTd  and  performing  a 
truncated  singular  value  decomposition  of  ZTZ.  The  level  of  truncation  was  chosen  such  that  singular  values 
less  than  the  floating  point  precision  of  MatLab  were  omitted.  Figure  3.4  contains  the  results  of  fitting  a  single 
channel  of  TEM  data  (Geonics  EM-63  sensor)  recorded  on  a  test  stand  over  a  40  mm  projectile  using  both  the 
patch  model  (Figure  3.4(a))  and  the  ring  model  (Figure  3.4(b)).  The  projectile  was  oriented  at  a  45  degrees  dip 
and  its  center  located  60  cm  below  the  sensor.  Both  models  produce  excellent  fit  to  the  data  (as  is  indicated  by 
correlation  coefficients  close  to  unity  and  small  misfit).  The  model  with  patches  fits  the  data  particularly  closely; 
misfits  are  five  orders  of  magnitude  less  than  the  data.  This  exceptional  fit  raises  the  issue  of  over-fitting  the  data 
by  fitting  the  noise.  Although  noise  levels  remain  low  within  the  well-controlled  environment  of  the  test  stand, 
a  close  look  at  the  low  left  panel  of  Figure  3.4(a)  suggests  the  presence  of  noise.  Indeed,  the  recorded  anomaly 
exhibits  more  shape  features  than  expected  for  a  simple -geometry  40  mm  projectile,  a  clear  effect  of  noise  adding 
high  frequency  variability  to  the  observation.  In  Figure  3.4(b),  discretizing  the  spheroid  surface  with  rings  of 
constant  charge  density  predicts  a  smooth  response  (low  center  panel);  this  might  be  considered  to  be  a  better 
solution  since  it  is  simpler  and  does  not  have  extreme  structure.  The  choices  of  rings  over  patches  and  the  number 
of  rings  are  first  steps  toward  a  regularization  of  the  surface  charge  model  imposed  by  discretization.  Shubitidze 
et  al.  (2005b)  restricted  the  number  of  rings  to  9.  We  tested  the  feasibility  and  effect  of  using  more  or  less  rings. 

In  Figure  3.4(b)  the  resulting  distribution  in  rings  remains  characterized  by  extreme  variations  between 
positive  and  negative  values  that  are  several  orders  of  magnitude  greater  than  the  magnitude  of  charge,  especially 
for  bottom  rings.  This  also  seems  to  create  additional  energy  outside  of  the  domain  used  for  the  inversion  (see 
higher  misfit  in  the  lowest  part  of  the  low  right  panel  of  Figure  3.4(b)).  Additional  tests  with  both  test  stand  and 
synthetic  data  show  that  solutions  for  charge  distribution  and  total  magnetic  charge  are  sensitive  to  variations  in 
target  positioning  and  orientation,  as  well  as  data  noise.  Higher  levels  of  regularization  are  therefore  warranted 
for  the  use  of  either  rings  or  patches. 
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(a)  Fitting  noisy  data  with  a  charge  density  distribution  on  the  surface  of  the  spheroid. 
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(b)  Fitting  the  same  noisy  data  as  in  (a)  with  a  ring  model.  The  ring  model  no  longer  over  fits  the  data  noise,  resulting  in  a  higher 
residual  error. 


Figure  3.4:  40  mm  projectile  (45  degrees  inclination,  60  cm  below  sensor).  Both  solutions  of  charge  distribution 
produce  large  oscillations. 
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Additional  characteristics  of  the  problem  include: 


1.  As  the  number  of  model  parameters  (rings,  patches,  Fourier  coefficients)  increases,  data  are  better  fit 
but  the  amplitude  variations  in  the  charge  density  of  adjacent  patches  or  rings  increases.  Also,  there  is  enhanced 
sensitivity  due  to  the  noise. 

2.  The  total  magnetic  charge,  the  correlation  coefficient  and  the  misfit  depend  on  the  size  and  shape  of 
the  spheroid  used  in  the  calculation.  For  instance,  large  elongated  items  like  a  105  mm  High  Explosive  Anti 
Tank  (HEAT)  round  would  be  better  represented  with  a  long  spheroid.  This  suggests  that  a  judicious  choice 
of  spheroid  should  be  made  before  carrying  on  the  inversion.  Practically,  this  would  imply  that  a  predefined 
and  limited  (3-5)  set  of  spheroids  be  chosen  to  test  possible  forward  models.  The  spheroid  should  not  be  too 
small  because  this  has  a  negative  effect  on  the  condition  number  of  Z:  physically,  a  smaller  spheroid  tries  to 
accommodate  all  the  features  of  the  target,  yielding  large  contrast  among  surface  charges. 


3.5  Properties  of  the  Normalized  Surface  Magnetic  Charge 

The  previous  section  showed  that  direct  inversion  of  EM  data  does  not  allow  recovering  of  a  stable  surface  charge 
distribution  and  that  a  regularization  of  inversion  is  needed.  In  order  to  guide  the  regularization  process  with 
physical  considerations,  we  now  re-examine  the  assumption  that  surface  magnetic  charges  can  be  normalized  to 
produce  a  surface  distribution  that  depends  on  the  object  only  and  is  independent  of  excitation  (relative  position 
and  orientation  of  object  and  illuminating  field). 

3.5.1  Field  on  a  surface  with  charges 

The  field  on  the  surface  of  a  charge  distribution  can  be  calculated  by  integrating  equation  3.2  by  using  a  Gaussian 
pillbox  on  a  charge  surface  (Figure  3.5).  The  integration  over  the  volume  gives 


(a)  Gaussian  pill  box  on  a  charge  distrib- 
tion  crm. 


(b)  Integrating  the  singularity  by  distorting 
the  charge  surface  with  a  semihemispheri- 
cal  surface. 


Figure  3.5:  Geometries  for  determining  the  field  on  a  charge  surface. 


(3.18) 


The  divergence  theorem  leads  to  surface  integrals,  and  the  following  expression 


(3.19) 


By  symmetry,  we  can  write 


(3.20) 


Hence  the  charge  distribution  is  proportional  to  the  normal  component  of  the  field  at  the  surface. 
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3.5.2  The  normalized  surface  magnetic  charge  on  a  spherical  surface  for  a  dipole  field 


To  understand  the  properties  of  the  normalized  surface  magnetic  charge,  it  is  instructive  to  investigate  the  nor¬ 
malized  charge  distribution  behavior  for  a  dipolar  field.  The  analytic  solution  for  the  secondary  field  of  a  sphere 
in  a  uniform  field  is  a  dipole.  Consider  a  permeable  and  conducting  sphere  of  radius  a  illuminated  by  a  uniform 
primary  field  Hp.  At  a  time  t  =  0  the  primary  field  is  terminated,  and  eddy  currents  are  induced  in  the  sphere; 
they  subsequently  decay  because  of  the  finite  conductivity  of  the  sphere.  The  secondary  field  Hsc  generated  by 
the  decaying  currents  is  dipolar: 

Hsc  (*)  =  4^3  m  (t)  ■  (3ff  -  i)  (3.21) 

where  m  (t)  is  the  dipole  moment  induced  at  the  center  of  the  sphere  at  time  t,  r  is  the  distance  between  the 
observation  point  and  the  sphere  center,  r  is  the  unit  vector  pointing  from  the  sphere  center  to  the  observation 
point  P,  and  I  is  the  identity  dyadic.  The  dipole  moment  is 

0-7 r 

m  (t )  =  —  Hp  L  ( t )  (3.22) 

Mo 


where 


OO 

L  ( t )  =  6a3 

S=1 


exp  (— f/r) 

Qs  +  {/J-r  —  1)  (Mr-  +  2) 


(3.23) 


with  r  =  a/ia2 /q2  (Kaufman  and  Keller,  1985).  The  values  qs  are  roots  to  the  transcendental  equation 


tangs 


(Mr  -  1)  qs 
Qs  +  (Mr-  -  1) ' 


(3.24) 


Equation  3.21  reveals  that  the  secondary  field  of  a  sphere  in  a  uniform  primary  field  is  equivalent  to  the  field  of  a 
single  magnetic  dipole  located  at  the  center  of  the  sphere  and  oriented  parallel  to  the  primary  field.  The  size  and 
material  properties  (i.e.  conductivity  and  magnetic  permeability)  are  contained  within  the  function  L  (t).  Now  let 
us  enclose  the  sphere  with  a  spherical  charge  surface  S  that  has  a  radius  ag  (Figure  3.6(a)).  On  the  surface  S,  the 


Figure  3.6:  Sphere  surface, 
secondary  field  produced  by  the  sphere  will  be 


H(f,re  S) 


2nag 


m  (t) 


- =-HpL  (f) 

Mo<4 


(3.25) 
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The  surface  charge  distribution  can  be  expressed  as  a  function  of  the  primary  field  by  taking  the  normal  of  H  from 
equation  3.25,  and  substituting  into  equation  3.20 


<xm  —  2/r0  H  •  n 

=  2  -U(t)  [Hp  •  n] 
as 


(3.26) 


The  NSMC  is  then 

qm=2-^L(t)  (3.27) 

aS 

and,  from  equation  3.6,  the  TNMC  is  the  surface  integral  of  qm  on  S 

Q=—L  (f)  (3.28) 

as 

Equation  3.28  shows  that,  for  the  case  of  a  sphere  in  a  uniform  field,  the  NSMC  is  uniform  on  the  surface  S. 
The  normalized  magnetic  charge  is  proportional  to  L  (t)  and,  therefore,  the  normalized  charge  for  a  sphere  is 
a  measure  of  the  size  and  material  properties  of  the  sphere.  However,  the  total  normalized  magnetic  charge  is 
inversely  proportional  to  as  and  is  thus  dependent  on  how  we  choose  our  charge  surface  S. 

The  above  calculations  demonstrate  that,  for  a  fixed  ag,  the  TNMC  is  a  reasonable  discriminant  for  modeling 
spheres.  This  is  due  to  the  geometry  of  the  problem  which  causes  the  induced  dipole  to  be  parallel  to  the  primary 
field.  Let  us  consider  a  non-spherical  target  whose  response  is  generated  by  charges  on  a  spherical  surface  (for 
example.  Figure  3.6(b)).  The  field  of  an  compact  metallic  target,  can  be  approximated  by  a  dipole  model.  We  can 
write  the  induced  dipole  as 


m  =  M  •  Hp 

where,  following  the  coordinate  system  of  Figure  3.6(b), 
27T 


M  = 


Mo 


Lx  (f)  0  0 

0  Ly  ( t )  0 

0  0  Lz  (t) 


(3.29) 


(3.30) 


For  a  sphere  L  (t)  =  Lx  (t)  =  Ly  (t)  =  Lz  ( t ).  The  charge  distribution  for  the  dipole  model  can  be  written  as 
<xm  —  2/r0  H  •  n 


[M  •  Hp]  •  n 


(3.31) 


7r  a' 


and  the  NSMC  is 


tio  [M  ■  Hp]  ■  n 
<?”1  7tac  Hp  •  n 


(3.32) 


It  is  clear  from  Equations  3.31  and  3.32  that  charge  distribution  <jm  and  normalized  charge  distribution  qrn 
is  dependent  on  the  direction  of  excitation  Hp  and  the  size  of  the  charge  surface.  For  example,  for  a  unit  vertical 
primary  field  Hp  =  z  the  normalized  charge  distribution  will  be  uniform: 


Qrri  —  3  L/z  ( t ) 


with  a  total  normalized  charge  of 

Q  =  —Lz  ( t ) 
as 


(3.33) 


(3.34) 
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If  the  primary  field  is  Hp  =  (k  +  y)  /y/2,  then  the  normalized  charge  distribution  will  vary  along  the  azimuthal 
direction  0  <  A  <  2ir: 


Qrri  (A)  ^ 

aS 

The  TNMC  is  then 

47 T 

Q=  —  [LX  ( t )  +  Ly  (i)  (3.36) 

as 


Lx  ( t )  cos  (A)  +  Ly  (t)  sin  (A) 


cos  (A)  +  sin  (A) 


(3.35) 


Clearly,  (1)  the  TNMC  is  dependent  on  the  radius  of  the  spherical  charge  surface,  and  (2)  the  TNMC  is 
dependent  on  the  direction  of  the  primary  field  that  illuminates  the  target. 


3.5.3  Investigating  the  normalized  magnetic  charge  using  the  Method  of  Moments 

In  order  to  study  the  normalized  magnetic  charge  on  a  spheroidal  surface,  we  will  solve  the  modeling  integral 
equations  with  the  MoM. 


Define  integral  equation 

Let’s  determine  the  field  on  the  surface  S.  There  is  a  singularity,  therefore  we  use  a  principle  value  integral 
and  explicitly  calculate  the  contribution  of  the  singularity. 


Hs  (r)  = 


lim 


1 


4t THo  Js_se  |r  -  r 


Tig  crm  (V)  d,S'  +  lim 


r  —  r 


+0  47 T[i0 


—  om  (r)  dS'  (3.37) 


where  Se  is  the  hemisphere  in  Figure  3.5(b).  The  singularity  as  r  — »  r'  is  integrable,  and  can  be  determined  by 
using  the  geometry  of  Figure  3.5(b)  and  letting  e  — >  0.  The  value  of  the  field  at  the  surface  is 

Hs  (r)  =  lim  — * —  f  4 — *-3  (r')  dS' 

e^o  47 THo  Js  |r  -  r'|3 


lim  — * —  f  ^  am(r')  dS' 
Air  Ho  JSe  e3 


2^n(r)  am(r) 


Therefore,  the  field  is 


(3.38) 


Hs  (r)  =  lim 


1 


L.w^r,[r')dS'+  L;n(r)  <,”(r) 


+0  AttHo 

where  r  £  S.  We  want  to  match  the  normal  components  of  the  magnetic  flux  B  at  the  surface: 


rP.V. 


n  (r)  •  (r  —  r') 


1 


r'|3 


Pm  (r')  dS'  L  -  Pm  (r) 


(3.39) 


(3.40) 


The  equations  for  the  normalized  surface  magnetic  charge  are  then 


rP.V. 


n  (r)  •  (r  —  r') 


r'|3 


qm  (r')  [Hp  (r')  •  n  (r')]  dS'+ 


-qm  (r)  [Hp  (r)  •  n  (r)]  (3.41) 
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Define  Basis  for  qm 


The  first  step  to  solving  for  qm  is  to  partition  the  surface  into  N  patches  of  area  A Sn.  We  assume  pulse 
basis  functions  to  represent  q.m: 


N 


Qm  (**)  —  ^  '  Grifn 


where 


fn  = 


Substitution  gives 


1  for  r  €  ASn 
0  otherwise 


(3.42) 


(3.43) 


„(r).B‘«T  ^  "(r)'(r-r') 

V  '  A  srr- 


47 r 


r  —  r 


•/ 1 3 


N 


Y  anfn  (r') 


n=l 


Hp  (r')  ■  n  (r')]  dS'+ 


N 


^  ^  Q-nfn  (* 


.n— 1 


Since  we  have  pulse  basis  functions,  we  can  approximate  the  integral  as 

Is  5  7s,, 


Hp  (r)  •  n  (r)]  (3.44) 


(3.45) 


n— 1  1 


which  gives 


n(r).B 'Af 

47rfrf4AS„  |r-r'|3 
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Y  ak^k  (r/) 


Lk—l 


Hp  (r')  •  n  (r')]  dS'+ 


N 


y '  cLnfn  (i 


.n—  1 


Hp  (r)  •  n  (r)]  (3.46) 


Note  that  we  changed  indices  for  the  pulse  basis,  in  the  the  first  term  on  the  right  hand  side.  Now,  fj.  (V)  will  act 
like  a  delta  function,  i.e.  the  integral  will  be  zero  unless  k  =  n.  Therefore,  we  rewrite  the  integral  as 

“  (r)  ‘  BS  ~  h  Y  /  ,  i./|3r)a"  [RP  (r')  •  n  (r')]  dS'+ 

n=rAS» 


1 


N 


y  \  Ctnfn  (l 


Hp  (r)  •  n  (r)]  (3.47) 


Defining  weighting  functions 

By  defining  a  basis  function  for  qm,  Equation  3.47  represents  an  approximation  to  the  linear  operator  defined 
in  3.41.  If  we  define  the  exact  integral  equation  as  L qm  =  g,  then  we  have  a  non  zero  residual 


JV 


N 


R  =  L qm  -  g  =  L  anfnj  -  g  =  a„L/nJ  -  g  (3.48) 

Weighting  functions  Wrn  are  used  to  make  the  residual  zero  at  a  finite  number  of  points,  i.e., 

f(Wm,  R)  (3.49) 
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(3.50) 


N 


]Ta„(Wm,  L)  =  (Wm,g) 


n=l 

This  defines  a  new  linear  system  to  solve  for  the  basis  coefficients  an 
G  a  =  h 

We  choose  point  matching  weighting  functions 

Wm  =  6  (r  —  rm )  where  m  =  1, N 

and  define  an  inner  product 


{f,h)=  f  f{v')h{r')dS' 

Js 


(3.51) 


(3.52) 


(3.53) 


Applying  to  the  first  term,  gives 
(Wm,  1st  term)  = 

[  S  (r  -  rm) 
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=—Y> 

47 r 


n=l 
JV  „p.y. 


dS 


(3.54) 


n— 1 


/  ASn 


[Hp  (rf)  ■  —  (r')]  d5- 


Applying  the  weighting  to  the  second  term,  gives 


r  h  N 

(Wm, 2nd  term)  =  6( r  -  rm)  -  V]  anfn  (r)  [Hp  (r)  •  n  (r)] 

^  2  „=i 


dS 


AT 


=  9  X!  (rm)  [HP  ‘  n 


(3.55) 


n—1 

=  [Hp  (rm)  •  n  (rm)] 

We  now  have  the  elements  of  G: 

1  fP  V-  n  (rm)  •  (rm  -  r')  , 


Gmn  —  A 

47T  , 


A5n 


r'|3 


[Hp  (r')  •  n  (!•')]  dS'  +  l-5mn  [Hp  (rm)  ■  n  (rm)]  (3.56) 


The  surface  integration  in  each  patch  AS  is  approximated  by  assuming  that  all  values  are  constant  within  each 
patch  A Sn 

Gmn  —  “7  j  [H  (^*n)  ’  n  (1  0-mn)  H“ 

|r  m  rn| 

which  can  be  rewritten  as 

1  n  (rm)  •  (rm  -  r„) 


-dm„[Hp(rm)-n(rm)]  (3.57) 


Gr 


4tt  \rm  —  rn|3 
\  [HP  (rm)  •  n(rm)] 


[Hp  (r„)  •  n  (r„)]  A Sn  if  m  ^  n 
if  m  =  n 


(3.58) 


Now  for  the  right  hand  side. 

hm  =  (Wm,n  (r)  •  Bs  (r))  =  [  6  [( r  -  rm)  n  (r)  •  Bs  (r)]  dS 

Js 

=  n(rm)  •  Bs  (rm) 


(3.59) 
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Examples  of  Method  of  Moments  modeling 

We  apply  the  MoM  modeling  to  the  problem  geometry  of  Figure  3.7.  The  field  of  the  target  is  represented 
by  a  dipole. 


Location  of 
predicted  and 


Figure  3.7:  Geometry  of  example.  We  match  boundary  conditions  on  a  surface  and  compare  the  data  predicted  by 
the  charge  distribution  on  the  surface  with  the  data  predicted  by  a  dipole  model. 


For  this  test  we  will  use  the  method  of  moments  to  determine  the  charge  distribution  am  and  the  normalized 
charge  distribution  qm  on  a  surface  that  encloses  the  dipole.  A  separate  matrix  equation  is  solved  for  om  and  qm, 
i.e.,  we  do  not  solve  for  the  qm  (or  <rm)  charge  distribution  and  multiply  (or  divide)  by  the  normal  component  of 
the  primary  field  to  obtain  am  (or  qm).  We  will  assume  that  the  primary  field  Hp  is  uniform.  The  field  predicted 
by  the  charges  obtained  by  method  of  moments  and  the  field  predicted  by  the  dipole  are  compared  on  a  surface  1 
m  above  the  dipole  location. 


Example  1:  Solution  on  a  spherical  surface  The  first  example  is  for  a  sphere  in  a  uniform  held.  As  was 
outlined  in  the  Section  3.5.2,  the  dipole  moment  is  due  to  a  sphere  and  is,  therefore,  parallel  to  the  primary  held. 
For  this  example,  we  will  try  and  use  a  spherical  charge  surface.  The  charge  surface  is  discretized  uniformly  in 
azimuth  (<f>)  and  dip  (??)  angles. 

Figure  3.8(a)  plots  the  recovered  un-normalized  (<rm)  and  normalized  (qm)  charge  distributions  in  the  case 
where  the  primary  held  is  vertical.  The  radius  of  the  charge  surface  is  0.15  m.  As  predicted  by  the  analysis  in 
the  previous  section,  the  normalized  charge  is  essentially  uniform.  The  un-normalized  charge  is  symmetric  about 
z  =  0,  and  is  consistent  with  a  vertical  magnetic  dipole.  The  held  predicted  by  both  the  normalized  and  un- 
normalized  distributions  are  compared  to  the  held  predicted  by  a  dipole  at  the  center  of  the  sphere.  The  forward 
modeled  charge  and  dipole  model  match  very  well,  indicating  that  the  charge  and  dipole  models  are  equivalent 
in  this  case.  Table  3.1  summarizes  our  results  for  illuminating  the  sphere  in  different  directions,  and  a  couple 
of  different  coarseness  levels  of  discretization.  The  total  unnormalized  magnetic  charge  is  zero  (to  numerical 
precision),  and  the  total  normalized  magnetic  charge  is  similar  in  each  case.  This  numerical  result  is  consistent 
with  the  analytic  result  presented  earlier:  the  total  normalized  magnetic  charge  is  stable  to  illumination  direction 
for  targets  with  spherical  symmetry. 

Figure  3.9  has  a  band  of  normalized  magnetic  charge  that  is  zero.  This  is  due  to  the  normal  component  of  the 
horizontally-directed  primary  field  being  zero  at  this  point.  The  corresponding  elements  in  the  MoM  modeling 
matrix  also  go  to  zero  at  this  point,  making  the  determination  of  the  normalized  magnetic  charge  at  those  points 
poorly  determined.  The  coarser  discretization  of  Figure  3.12  does  not  have  surface  patches  that  have  normal 
perpendicular  to  the  primary  field.  Consequently,  all  the  surface  patches  are  well  determined. 
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Primary 

Field 

num 

</> 

num 

V 

Li 

L2 

total  unnormalized 
magnetic  charge 

total  normalized 
magnetic  charge 

Figure 

Number 

Vertical 

30 

30 

1 

1 

-4.3e-014 

20.6 

3.8 

Horizontal 

30 

30 

1 

1 

-6.2e-015 

19.3 

3.9 

Vertical 

20 

20 

1 

1 

-1.4e-013 

20.91 

3.10 

45  degrees 

20 

20 

1 

1 

-9.7e-014 

20.88 

3.11 

Horizontal 

20 

20 

1 

1 

-4.2e-016 

20.88 

3.12 

Table  3.1: 

Magnetic  charge  results  for  the  sphere  tests. 

Primary 

num 

num 

total  unnormalized 

total  normalized 

Figure 

Field 

<t> 

V 

Li 

L2 

magnetic  charge 

magnetic  charge 

Number 

Vertical 

21 

21 

4 

1 

-8.9e-014 

364.07 

3.13 

45  degrees 

21 

21 

4 

1 

-6.1e-014 

1052.86 

3.14 

Horizontal 

21 

21 

4 

1 

-9.1e-015 

22.00 

3.15 

Table  3.2:  Magnetic  charge  results  for  the  spheroid  tests. 


Example  2:  Solution  on  a  spheroidal  surface  Our  second  example  is  to  determine  the  charge  distribution  for 
an  axi-symmetric  target.  The  secondary  field  for  the  axi-symmetric  target  is  modeled  using  a  dipole  field  with 
L\  =  4  and  /,2  :  T.a  =  1.  For  this  example,  we  use  a  spheroidal  charge  surface.  The  spheroidal  surface  has  a 
length  of  40  cm  and  a  width  of  20  cm.  The  spheroidal  surface  is  oriented  with  the  major  axis  parallel  to  the  z-axis 
for  each  of  the  examples  of  this  section. 

Figure  3.13(a)  plots  the  recovered  un-normalized  ( <jm )  and  normalized  ( qm )  charge  distributions  in  the  case 
where  the  primary  field  is  vertical.  The  un-normalized  charge  is  symmetric  about  z  —  0,  and  is  consistent  with  a 
vertical  magnetic  dipole.  The  field  predicted  by  both  the  normalized  and  un-normalized  distributions  are  compared 
to  the  field  predicted  by  the  dipole  at  the  center  of  the  spheroid.  The  forward  modeled  charge  and  dipole  model 
match  very  well,  indicating  that  the  charge  and  dipole  models  are  equivalent  in  this  case. 

Table  3.2  summarizes  our  results  for  illuminating  the  spheroid  in  different  directions.  As  was  the  case  with 
the  sphere  examples,  the  total  un-normalized  magnetic  charge  is  zero.  In  each  case  the  charge  distribution  and 
dipole  fields  match  on  a  plane  1  m  away  from  the  center  of  the  spheroid.  In  this  case,  the  total  normalized  magnetic 
charge  is  not  independent  of  the  primary  field  direction. 
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(a)  Charge  Distribution  on  a  sphere.  The  top  panel  is  the  normalized  surface  magnetic 
charge  distribution  qm .  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge 
distribution  crm . 


(b)  Data  modeled  1  m  above  the  center  of  the  sphere 


Figure  3.8:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  sphere  with  a  vertical 
primary  field  excitation.  The  surface  was  discretized  with  900  patches. 
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(a)  Charge  distribution.  The  top  panel  is  the  normalized  surface  magnetic  charge  distri¬ 
bution  qm .  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge  distribution 
crm.  The  blue  strips  in  the  normalized  charge  distribution  correspond  to  points  where 
HP  •  n  =  0 


(b) 


Figure  3.9:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  sphere  with  a  hori¬ 
zontal  primary  field  excitation.  The  surface  was  discretized  with  900  patches. 
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(a)  Charge  Distributions.  The  top  panel  is  the  normalized  surface  magnetic  charge 
distribution  <?m.  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge  dis¬ 
tribution  crm. 


(b)  Data  modeled  1  m  above  the  center  of  the  sphere 


Figure  3.10:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  sphere  with  a  vertical 
primary  field  excitation.  The  surface  was  discretized  with  400  patches. 
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Normalized  charge 


(a)  Charge  distribution.  The  top  panel  is  the  normalized  surface  magnetic  charge  distri¬ 
bution  q-m .  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge  distribution 

&  m  ■ 
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(b)  Data  modeled  1  m  above  the  center  of  the  sphere 


Figure  3.11:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  sphere  with  a  primary 
field  excitation  at  a  45  degree  angle  (Hp  =  (y  +  z)  jyf%.  The  surface  was  discretized  with  400  patches. 
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(a)  Charge  distribution.  The  top  panel  is  the  normalized  surface  magnetic  charge  distri¬ 
bution  q-m .  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge  distribution 
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(b)  Data  modeled  1  m  above  the  center  of  the  sphere 


Figure  3.12:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  sphere  with  a  hori¬ 
zontal  primary  field  excitation.  The  surface  was  discretized  with  400  patches. 
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Normalized  charge 


(a)  Charge  distribution.  The  top  panel  is  the  normalized  surface  magnetic  charge  distri¬ 
bution  qm  •  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge  distribution 
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(b)  Data  modeled  1  m  above  the  center  of  the  spheroid 


Figure  3.13:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  spheroid  with  a 
vertical  primary  field  excitation. 
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Normalized  charge 


(a)  Charge  distribution.  The  top  panel  is  the  normalized  surface  magnetic  charge  distri¬ 
bution  q-m .  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge  distribution 
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(b)  Data  modeled  1  m  above  the  center  of  the  spheroid 


Figure  3.14:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  spheroid  with  a 
primary  field  at  a  45  degree  angle  (Hp  =  (y  +  z)  / -\/2).. 
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(a)  Charge  distribution.  The  top  panel  is  the  normalized  surface  magnetic  charge  distri¬ 
bution  q-m .  The  bottom  panel  is  the  un-normalized  surface  magnetic  charge  distribution 

&  m  ■ 


(b)  Data  modeled  1  m  above  the  center  of  the  spheroid 


Figure  3.15:  Method  of  Moments  solution  for  the  surface  magnetic  charge  distributions  on  a  spheroid  with  a 
horizontal  primary  field  excitation. 
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3.5.4  Comparison  to  SEA  modeled  secondary  fields 


The  previous  section  presented  the  distributions  of  surface  charges 


am  =  2p,0Hsc  •  n 


(3.60) 


and  normalized  surface  charges 


Hs 


Qm  —  ftQ 


HP  n 


(3.61) 


in  the  simplified  case  where  the  scattered  field  Hsc  could  be  represented  by  a  dipole  model  and  the  primary  field 
Hp  was  uniform.  Realistic  conditions  were  tested  for  a  cylinder  for  which  the  SEA  library  computed  scattered 
fields  at  the  surface  of  a  spheroid  enclosing  the  cylinder  and  with  dimensions  defined  by  the  user.  SEA  was  used 
after  verification  that  the  computation-intensive  MAS  and  simplified  SEA  codes  produced  the  exact  same  result. 
Realistic  primary  field  was  computed  by  modeling  the  effect  of  a  frequency  domain  square  transmitter  of  1  x  1  m. 


In  the  following,  only  the  real  part  of  the  scattered  field  measured  at  the  first  frequency  (0.1  mHz)  is  dis¬ 
played,  as  later  frequencies  and  the  imaginary  part  exhibit  similar  trends.  Each  figure  presents  on  a  spheroid 
encapsulating  the  same  cylinder  the  distribution  of  normal  component  of  the  scattered  and  primary  fields,  as  well 
as  the  normalized  charge  taken  as  their  direct  ratio  and  that  obtained  by  resolution  through  the  method  of  moment. 
The  effect  of  the  sizes  and  shapes  of  the  modeling  spheroid  (and  sphere)  and  those  of  the  relative  orientations  and 
positions  of  transmitter  and  target  were  tested. 


Reference  test 

As  a  reference  test  presented  in  Figure  3.16  the  cylinder  is  placed  vertically  directly  below  the  transmitter 
at  a  42.7  cm  depth.  The  enclosing  modeling  spheroid  has  semi-axes  of  length  21  cm  and  10  cm  and  a  resolution 
of  40  points  in  latitude  and  39  in  longitude,  with  longitude  ranging  from  0  to  27t  and  latitude  from  0  (North)  to  n 
(South).  In  each  panel,  quantities  are  shown  on  the  spheroid  by  unwrapping  its  surface  and  using  contours  with  a 
special  contour  in  magenta  to  indicate  zero  crossings  for  normalized  charge  density.  The  vertical  position  exactly 
below  the  transmitter  implies  that  all  fields  are  longitudinally  symmetric,  contour  lines  therefore  follow  latitudes. 
In  the  two  left  panels  of  Figure  3.16  the  normal  component  of  the  scattered  and  primary  fields  approach  zero  at 
different  latitudes  because  the  primary  field  is  not  uniform,  their  direct  ratio  therefore  yields  large  values  (infinite 
value  if  discretization  allowed)  in  the  upper  central  panel  while  MoM  gives  large  but  well  defined  values  (right 
panels).  Normalized  charges  computed  by  direct  ratio  or  MoM  show  similar  but  different  values,  their  integral 
over  the  surface  of  the  spheroid,  the  TNMC  differs  by  30%. 


Changes  of  orientation 

The  same  cylinder  and  its  spheroid  are  rotated  into  a  horizontal  position  in  Figure  3.17.  Although  distri¬ 
butions  of  normalized  charges  by  field  ratios  and  MoM  appear  qualitatively  similar,  with  most  of  the  spheroid 
covered  with  negative  charges,  differences  in  charge  amplitude  are  such  that  their  respective  Total  Magnetic 
Charge  (TMC)  obtained  by  numerical  integration  vary  by  a  factor  of  two  (with  the  resolution  applied  here).  Com¬ 
parison  with  the  previous  figure  shows  that  rotation  of  the  illuminating  field  by  90  °  stimulates  a  totally  different 
charge  distribution  at  the  surface  of  the  spheroid  and  TMC.  Dependency  on  the  illuminating  field  is  further  illus¬ 
trated  in  Figure  3.18  with  a  rotation  of  45  °  and  yet  another  normalized  charge  distribution  characterized  by  large 
regions  of  positive  and  negative  charges  and  a  different  value  for  TMC. 


Change  of  position 

The  transmitter  was  modeled  at  30  cm  to  the  left  of  the  center  of  the  cylinder  and  spheroid  for  the  tests 
presented  in  Figures  3.19-3.21.  For  each  case  the  normalized  charge  distribution  and  TMC,  derived  by  direct 
ratio  and  MoM,  show  different  patterns  and  values  and  confirm  strong  dependency  on  the  relative  position  and 
orientation  of  transmitter  and  modeled  target. 
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Figure  3.16:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  spheroid  enclosing  a  vertical  cylinder.  Reference  case. 
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Figure  3.17:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  spheroid  enclosing  a  horizontal  cylinder. 
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Figure  3.18:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  spheroid  enclosing  a  cylinder  at  45  °. 
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Figure  3.19:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  spheroid  enclosing  a  vertical  cylinder  located  30  cm  off  the  center  of  the  transmitter. 
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Figure  3.20:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  spheroid  enclosing  a  horizontal  cylinder  located  30  cm  off  the  center  of  the  transmitter. 
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Figure  3.21 :  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  spheroid  enclosing  an  oblique  cylinder  located  30  cm  off  the  center  of  the  transmitter. 


61 


Effect  of  resolution,  size  and  shape 


In  order  to  check  whether  the  observed  trends  are  artifacts  of  the  numerical  implementation  of  the  problem, 
effect  of  resolution  in  discretization,  size  and  shape  of  the  spheroid  were  tested  in  Figures  3.22-3.24.  In  all  cases 
the  charge  distribution  remains  qualitatively  similar.  The  total  magnetic  charge,  however,  varies  by  a  factor  60 
in  the  first  case  with  lower  resolution,  by  factor  -3  in  the  second  case  with  semi-axes  increased  by  50%  and  by  a 
factor  6  in  the  third  case  of  a  spherical  surface. 
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Figure  3.22:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  spheroid  enclosing  a  vertical  cylinder  using  lower  resolution. 


Result 

The  normalized  surface  charge  density  distribution  varies  at  the  surface  of  a  sphere  or  spheroid.  Its  sign 
varies  in  a  predictable  manner  that  reflects  the  projections  of  the  primary  and  secondary  fields  at  the  surface  of  the 
spheroid.  Moreover,  singularities  in  the  normalized  charge  distribution  can  occur  when  the  normal  component  of 
the  primary  field  does  to  zero.  The  singularities  are  avoided  only  if  the  normal  component  of  the  secondary  field 
also  vanishes  at  the  exact  same  location.  That  distribution  changes  when  the  relative  position  of  the  spheroid  and 
transmitter  varies,  therefore  the  normalized  surface  charge  density  distribution  is  not  a  characteristic  of  an  object. 
TMC  changes  with  the  position  and  orientation  of  the  object,  as  well  as  with  the  shape,  size  and  resolution  used 
for  the  spheroid,  an  effect  of  the  aforementioned  singularities.  This  result  is  in  direct  agreement  with  the  previous 
section,  this  time  using  spatially  varying  primary  field  and  a  physical  model  of  the  secondary  field  derived  from 
MAS. 


3.5.5  Conclusion 


The  series  of  tests  presented  in  this  section  has  shown  that  the  normalized  magnetic  charge  qm  is  a  function 
of  the  geometry  of  the  measurement.  In  order  for  the  total  normalized  magnetic  charge  Q  to  be  an  effective 
discriminant,  we  need  to  ensure  high  quality  data  is  taken  where  the  target  is  excited  at  numerous  angles  of 
primary  field  excitation.  The  total  magnetic  charge  for  such  a  measurement  then  reflects  an  “average”  of  the  total 
magnetic  charges. 
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Figure  3.23:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  large  elongated  spheroid  enclosing  a  vertical  cylinder. 
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Figure  3.24:  Distribution  of  the  normal  component  of  the  scattered  and  primary  fields  and  their  ratio  at  the  surface 
of  a  sphere  enclosing  a  vertical  cylinder. 
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3.6  Regularized  Inversion  for  a  Surface  Charge  Distribution 


Section  3.3.1  describes  the  numerical  forward  model  that  calculates  the  secondary  field  produced  by  a  surface 
charge  distribution  expanded  by  a  pulse  basis  function.  For  a  given  charge  distribution  surface,  the  forward  model 
is  linear  and  can  be  written  as  Zm  =  d.  The  objective  of  the  inverse  problem  is  to  determine  a  charge  distribution 
m  from  a  data  set  d.  There  are  a  number  of  challenges  in  determining  an  optimal  m.  First,  the  number  of 
parameters  used  to  represent  the  charge  distribution  (i.e.,  the  coefficients  multiplying  the  basis  functions)  can  be 
greater  than  the  number  of  data.  In  such  a  case  the  problem  is  under-determined.  Second,  data  are  noisy  and 
special  care  has  to  be  taken  to  ensure  that  our  model  does  not  fit  the  noise.  Third,  the  data  may  be  insufficient  to 
constrain  all  the  model  parameters.  For  example,  if  the  illuminating  field  is  only  in  a  single  direction,  then  the 
problem  is  ill-posed. 

As  demonstrated  in  Section  3.4,  determining  the  charge  distribution  by  minimizing  a  data  misfit  objective 
function  is  ill-advised.  Figures  3.4(a)  and  (b)  demonstrated  that  minimizing  the  data  misfit  ||Zm  —  d|j  produces 
large  spikes  in  the  charge  distribution  which,  in  turn,  lead  to  very  unstable  estimates  of  the  total  normalized  mag¬ 
netic  charge.  In  order  to  solve  such  an  ill-posed  inverse  problem,  the  data  fitting  problem  has  to  be  reformulated 
in  order  to  introduce  prior  information  of  the  model  that  we  seek.  There  are  several  ways  to  do  this  reformulation. 
A  statistical  framework  (Tarantola,  1987;  Scales  and  Tenorio,  2001)  is  appropriate  when  our  notion  of  the  charge 
distribution’s  characteristics  can  be  represented  by  probabilities.  The  statistical  approach  can  be  used  to  determine 
a  single  model  by  determining  the  maximum  a  posteriori  probability  model,  which  can  be  determined  by  solving 
an  optimization  problem. 

In  the  previous  section  we  saw  that  the  surface  integral  of  the  un-normalized  surface  charge  (am)  was  zero. 
For  the  remainder  of  this  report,  we  will  use  the  terms  TMC  and  TNMC  to  represent  the  same  quantity,  i.e.,  Q. 


3.6.1  Formulation  of  the  regularized  inverse  problem 


We  use  an  alternative  approach  where  the  prior  information  is  introduced  through  the  introduction  of  a  regulariza¬ 
tion  functional,  or  model  objective  function  (Parker,  1994;  Menke,  1989).  The  regularization  functional  quantifies 
model  features  (such  as  amplitude,  smoothness,  and  curvature),  and  provides  a  quantitative  means  of  choosing 
models  that  are  consistent  with  our  preconceived  notion  of  the  characteristics  of  the  model.  We  achieve  this  by 
solving  an  optimization  problem: 

minimize  $  =  $d  +  A$m  (3.62) 

where  <l\i  is  the  data  misfit  objective  function,  <f>m  is  the  model  objective  function,  and  0  <  A  <  oo  is  the 
regularization  parameter.  The  data  misfit  function  measures  the  misfit  between  data  predicted  by  a  model  and  the 
observed  data  that  we  seek  to  fit.  The  relationship  between  the  modeled  and  observed  data  can  be  written  as 

dobs  =  Zm  +  noise  (3.63) 

Our  goal  is  to  fit  the  data  without  fitting  ’’noise”,  where  we  define  ’’noise”  as  any  component  of  the  observed  data 
that  should  not  be  accounted  for  by  the  modeling.  Therefore  noise  would  include  sensor  measurement  errors, 
and  modeling  errors.  For  our  work  we  will  assume  that  the  "noise”  is  Gaussian  and  independent.  With  this 
assumption,  we  can  write  the  data  misfit  as 

$d  =  ||Wd(Zm-d)||i  (3.64) 

where  Wd  is  the  diagonal  matrix  Wd  =  diag  (1/ei)  / y/(lV),  with  e,  defined  as  the  standard  deviation  of  the  ith 
datum  and  N  the  number  of  data  points.  With  this  definition  of  W,/  the  expectancy  of  data  misfit  /:’[<k|]  =  1  if 
the  model  is  consistent  with  the  standard  deviation  of  data. 

The  model  misfit  $m  is  defined  as 

$m=||Wm(m-m0)||2  (3.65) 
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where  Wm  is  the  regularization  operator  and  m0  a  reference  model. 


Having  observed  in  Section  3.5  that  for  a  homogeneous  sphere  the  normalized  charge  is  uniform,  we  assume 
at  first  order  that  the  charge  distribution  would  be  uniform  (i.e.,  similar  to  a  sphere),  and  that  additional  structure 
to  the  charge  distribution  would  be  allowed  in  order  to  fit  the  data.  Therefore  we  define  the  reference  model  m0 
to  be  the  uniform  charge  distribution,  where  the  charge  density  on  each  surface  element  of  the  discretization  can 
be  defined  as: 


[m0]  =  qunl  =  median 


N 


h-vE^ 

j= i 


(3.66) 


Because  un-regularized  inversion  of  data  such  as  in  Figure  3.4(b)  shows  large  oscillations  of  charge  amplitude 
between  neighboring  rings,  we  choose  to  impose  a  smoother  model  and  define  Wm  as  the  first  derivative  operator. 
We  note  that  in  that  case  Wmm0  =  0;  m0  remains  useful  to  define  an  approximation  to  the  total  magnetic  charge 
of  the  uniform  distribution,  while  the  true  uniform  distribution  is  identified  after  solving  the  optimization  problem 
of  Equation  3.62  and  taking  the  most  regularized  model,  the  smoothest  model  being  a  uniform  one.  This  subtlety 
arises  from  our  finding  that  these  two  definitions  of  the  uniform  model  can  depart  from  one  another  when  noise 
levels  are  high,  thus  introducing  a  bias  in  Equation  3.66  through  the  data  Hsc. 


The  parameter  A  is  chosen  to  balance  the  trade-off  between  the  fit  to  the  data  and  the  a  priori  information 
introduced  through  the  regularization.  There  are  several  methods  for  determining  the  trade-off  parameter  (see 
Farquharson  and  Oldenburg,  2004,  and  references  therein).  For  this  work  we  adopt  the  L-curve  criterion  for 
determining  the  trade-off  parameter.  The  L-curve  technique  involves  generating  an  L  —  curve  or  Tikhonov  curve 
(<I>d  =  /  (4*111)).  which  are  the  values  of  4\|  and  <f>m  for  different  values  of  the  trade-off  parameter.  The  optimal 
A  is  chosen  to  be  at  the  corner  of  the  ”L”  in  the  Tikhonov  curve,  i.e.,  the  point  with  maximum  curvature.  By 
adopting  this  criterion,  the  data  misfit  is  large  enough  to  accommodate  the  presence  of  noise,  while  being  small 
enough  to  provide  a  satisfactory  solution  to  Zm  =  d. 


We  can  now  write  the  inverse  problem  as 

minimize  <f>  =  ||W„  (Zm  -  d)  |||  +  A||Wm  (m  -  m0)  ||| 

The  solution  is  obtained  by  taking  the  gradient  of  the  $  and  setting  V<f>  =  0: 

mrec  =  (ZTWjWjZ  +  AW^Wm)_1  (ZTWjWjZd  +  AW^Wmm0 
where  we  define  mrec  to  be  the  recovered  model. 


(3.67) 

(3.68) 


Figure  3.25  demonstrates  the  effect  of  a  regularization  on  the  ring  model  with  the  same  noisy  data  as  in 
Figure  3.4(b).  The  regularized  charge  model  shows  a  smooth  distribution  of  charge  density  along  the  spheroid,  in 
sharp  contrast  with  the  large  amplitudes  observed  in  the  un-regularized  case.  Additional  experiments  with  20,  40 
and  60  rings  show  that  the  same  regularization  procedure  yields  the  exact  same  distribution.  For  the  remainder 
of  this  study  we  use  40  rings  to  discretize  the  spheroidal  surface,  as  we  feel  that  this  level  of  discretization 
is  sufficiently  fine  to  model  the  charge  distributions  of  relatively  complex  objects  without  requiring  excessive 
computational  effort. 

The  low  data  misfit  and  high  correlation  coefficient  between  the  observed  data  and  the  model  prediction 
are  similar  to  those  of  the  un-regularized  model.  However,  the  calculated  total  magnetic  charge  derived  from  the 
different  models  differ  by  several  orders  of  magnitude.  The  patch  model,  the  unregularized  ring  model  and  the 
regularized  ring  model  have  total  magnetic  charge  values  of  38  x  103,  — 1.09  x  101  and  —2.15  x  10-2,  respectively. 
For  the  total  magnetic  charge  to  be  a  stable  discrimination  criterion  we  must  understand  how  discretization  and 
regularization  affect  the  total  magnetic  charge. 

Figure  3.26  demonstrates  the  effect  of  regularization  on  the  amplitude  of  the  charge  density  distribution. 
At  early  stages  of  the  regularization  (small  A),  charges  vary  between  104  and  — 104  (the  top  panel).  At  later 
stages,  larger  values  of  A  bring  down  the  amplitude  by  a  factor  over  105,  shift  the  peaks  of  oscillations  and  force 
convergence  toward  a  flat  uniform  distribution  near  the  value  of  m0. 
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Figure  3.25:  Regularized  ring  solution  for  40  mm  projectile  (45  degrees,  60  cm  below  sensor).  In  top  central 
panel,  crosses  show  the  10-ring  solution,  the  solid  red  line  the  40-ring  solution  using  the  exact  same  automated 
regularization  procedure  (described  in  section  3.6.3).  Misfit,  expectancy,  correlation  and  total  magnetic  charge 
exactly  coincide  for  10-40  rings. 


3.6.2  Total  magnetic  charge  and  regularization 


The  left  panel  of  Figure  3.27  illustrates  the  concept  of  L-curve  for  the  previously  exposed  case  of  a  40  mm 
projectile  at  45  degrees  and  60  cm  below  an  EM-63  sensor.  A  common  method  for  choosing  an  appropriate 
regularization  is  to  pick  the  one  that  corresponds  to  the  corner  of  the  L-curve.  As  regularization  becomes  larger 
the  misfit  grows  and  the  charge  distribution  becomes  smoother.  Conversely,  for  weak  regularization  the  model 
can  fit  the  data  so  well  that  it  also  fits  noise,  and  therefore  becomes  sensitive  to  noise.  The  corner  of  the  L-curve  is 
shown  with  a  star.  The  large  dot  represents  an  alternative  pick  explained  in  the  next  section.  The  right  panel  shows 
L-curves  for  the  same  projectile  at  different  depths  and  orientations.  For  instance,  the  dashed  line  corresponds 
to  a  vertical  orientation,  where  regularization  plays  an  important  role  on  model  misfit,  as  reflected  by  the  larger 
variation  of  <I>d- 

The  previous  sections  showed  that  the  total  magnetic  charge  can  depend  on  discretization  and  regularization. 
This  first  result  is  troubling  because  it  contradicts  the  claim  of  Shubitidze  that  the  total  magnetic  charge  can  be 
used  as  a  classification  criterion  among  UXO,  and  therefore  that  it  should  remain  independent  of  object  position, 
orientation  and  background  noise.  Focusing  on  the  method  of  rings,  we  study  in  Figure  3.28  the  total  magnetic 
charge  Q  during  regularization  (parameter  A)  for  the  same  40  mm  target  as  in  Figure  3.27.  In  the  left  panel  the 
object  is  at  45  degrees.  Its  total  charge,  the  blue  solid  dotted  line,  varies  significantly  as  A  increases  and  seems 
to  converge  at  a  later  stage  toward  the  horizontal  dashed  line,  the  total  magnetic  charge  Quni  for  the  model  with 
uniform  distribution  of  charge  density.  This  is  expected  because  the  limit  case  of  regularization  by  smoothing  is 
a  flat  uniform  model.  As  in  Figure  3.27,  the  star  indicates  the  model  chosen  for  the  corner  of  the  L-curve  (Areg), 
the  dot  an  alternative  regularization.  Neither  of  these  points  marks  any  particular  position  for  the  Q  curve. 

We  consider  different  orientations  and  depths  of  the  same  target  in  the  right  panel  of  Figure  3.28  to  further 
study  the  total  magnetic  charge  of  the  40  mm  projectile.  Each  line  style  corresponds  to  the  same  cases  as  in 
Figure  3.27.  Figure  3.28  shows  that  each  configuration  of  depth  and  orientation  generates  different  variations  of 
Q  during  the  regularization  process.  These  variations  are  larger  when  the  inclination  of  the  object  is  greater,  as 
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Figure  3.26:  Charge  distribution  on  40  rings  along  modeling  spheroid  for  same  40  mm  item  (Latitude  +90  for 
the  pole  closest  to  the  surface,  0  equator,  -90  bottom).  Top  panel:  all  regularized  models,  least  regularized  with 
largest  amplitude,  A  in  log  scale  from  10_10  to  105.  Lower  panel:  sampled  models  for  lambda  =  4,  24,  123,  641, 
3350,  17500  and  above.  The  thick  dashed  line  corresponds  to  the  regularization  applied  for  Figure  3.25. 


Figure  3.27:  Example  of  L-curves  (i.e.  <I\|  vs.  <l>ln  curves)  for  a  40  mm  projectile.  Left  panel:  object  at  45 
degrees,  60  cm  below  sensor.  Stars  indicate  maximum  curvature,  large  dots  alternative  pick  for  regularization. 
Right  panel:  black,  red  and  green  curves  correspond  to  shallow  depth  (20-25  cm),  the  remaining  ones  to  deeper 
position  (60  cm);  black  solid  line  and  blue  line-dot  at  45  degrees,  red  dashed  and  cyan  dash-dot  vertical,  green 
crosses  and  magenta  plus  markers  horizontal. 
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illustrated  by  the  red  and  cyan  curves  that  span  several  orders  of  magnitude  with  change  in  sign,  as  opposed  to  the 
more  contained  variations  for  a  horizontal  position  (green  crosses  and  magenta  plus  signs).  Furthermore,  none  of 
these  Q  curves  seems  to  share  any  common  global  or  local  extremum,  which  a  priori  complicates  the  choice  of 
Q.  Star  markers  show  that  choices  of  Q(  Areg)  from  the  corner  of  the  L-curve  do  not  provide  a  consistent  value, 
thus  a  traditional  approach  to  regularization  is  not  compatible  with  a  unique  and  invariant  total  magnetic  charge 
of  the  40  mm  projectile. 

Searching  for  common  features  among  Q  curves,  we  consider  the  total  magnetic  charge  Q un;  for  the  model 
of  uniform  distribution  of  charge  density.  The  dashed  horizontal  lines  in  the  right  panel  of  Figure  3.28  show 
that  its  value  remains  similar  for  all  these  different  positions  and  orientations,  with  mean  -0.0231  and  standard 
deviation  0.0024.  Could  the  uniform  distribution  provide  an  adequate  forward  model  to  predict  data  recorded  at 
the  surface?  Detailed  analysis  reveal  that  these  models  are  limited  to  predicting  radial  EMI  responses,  yielding 
misfits  that  are  generally  too  large  to  qualify  them  as  acceptable  solutions  to  Zm  =  d.  This  is  particularly 
obvious  for  heterogeneous  objects  that  lie  in  sub-horizontal  position  because  their  EMI  response  reflects  their 
spatial  heterogeneity.  This  test  and  many  more  with  real  data  (some  of  which  are  discussed  in  the  next  section) 
and  synthetic  data  show  that  the  uniform  distribution  is  a  stable  characteristic  of  an  item  but  a  poor  forward  model. 
An  alternative  choice  of  regularization  is  thus  warranted  for  the  forward  model. 


Figure  3.28:  Total  magnetic  charge  vs.  regularization  parameter,  same  legend  and  setting  as  previous  figure.  Left 
panel  shows  large  variations  of  total  magnetic  charge  Qsmr  during  the  regularization.  The  dashed  horizontal  line 
shows  Qsinc(riio)  (also  called  Quni)-  Right  panel  shows  Qsmc  and  QUni  for  different  orientations  and  depths.  Quni 
is  similar  for  all  cases. 


3.6.3  Regularizing  for  a  stable  total  magnetic  charge 


The  previous  section  has  shown  that  the  corner  of  the  L-curve  does  not  guarantee  an  invariant  total  magnetic 
charge  whereas  a  uniform  charge  distribution  provides  a  stable  Q  but  a  potentially  deficient  forward  model,  espe¬ 
cially  for  sub-horizontal  items.  In  that  context,  we  consider  the  best  solution  to  be  one  that  satisfies  both  Q(A)  is 
close-enough  to  Quni  and  the  misfit  is  of  the  same  order  as  the  noise.  Theory  and  extensive  tests  with  real  mea¬ 
surements  and  synthetic  data  show  that  there  is  always  a  solution  to  the  first  condition  because  large  regularization 
imposes  smoothness  and  eventually  flatness  to  the  charge  distribution.  Besides,  there  are  often  models  whose  total 
magnetic  charge  approach  Qunl  at  early  stages  of  regularization,  as  illustrated  in  Figure  3.28  by  crossings  of  Qlnii 
curves  by  Q  curves  for  small  values  of  A.  These  models  are  potential  candidates  because  their  misfit  is  low  without 
them  fitting  noise  owing  to  some  regularization.  Other  possible  candidates  appear  as  regularization  increases  and 
Q  gradually  approaches  Qllni,  before  the  charge  distribution  becomes  uniform.  Practically,  quantitative  values 
must  be  assigned  to  properly  define  the  notion  of  approaching  C)m,j. 
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Method  1 


We  propose  to  choose  the  regularization  parameter  A  such  that 


||Q  (A)  —  Quni  II 

110  uni  || 


<  (5 


(3.69) 


and 

<  $d  <  K2  (3.70) 

Ai 

with  E  ($d)  =  1  if  the  model  fits  data  but  not  noise.  We  applied  this  method  by  first  setting  delta  =  20%,  K\  =  2 
and  K2  =  1.5,  allowing  an  error  of  20  %  on  Q  and  50%  on  noise  estimate.  If  there  is  no  solution  we  first  relax  6 
to  30%  and  set  K\  =  4,  K  >  =  2,  then  relax  further  to  <5  =  50%,  K\  =  8  and  K2  =  3. 

Method  2 

Noise  and  errors  from  instruments,  measurements  and  modeling  are  not  necessarily  straightforward  to  assess. 
A  priori  error  estimate  can  be  inaccurate  and  misleading  and  cause  the  method  above  to  fail,  even  though  we  find 
it  to  perform  well  most  of  the  time.  When  it  fails  we  can  either  choose  the  model  with  Q  closest  to  (A,lnj,  or  pick 
the  corner  of  the  L-curve.  In  either  case  we  loose  control  on  one  of  the  criteria.  To  circumvent  this  difficulty  one 
can  assume  that  the  corner  of  the  L-curve  provides  an  acceptable  misfit  thanks  to  its  balance  position  in  fitting  data 
and  noise,  with  expectancy  of  <f>d  is  f^($d)<I>d(Acorner)-  This  is  the  classical  justification  for  using  the  L-curve 
(e.g  Hansen,  1997;  Farquharson  and  Oldenburg,  2004).  Identifying  the  corner  of  the  L-curve  by  its  maximum 
curvature,  we  search  for  a  model  near  that  corner  and  with  a  total  magnetic  charge  close  to  Quni-  Similarly  with 
the  previous  method,  our  search  algorithm  identifies  the  regularization  parameter  A  that  simultaneously  minimizes 
the  relative  difference  between  Qun;  and  Q  on  the  one  hand,  and  between  L’LLi  j  and  d>d  on  the  other  hand: 

DQ  =  \max(Q(X),Quni)/min(Q(\),Quni))  —  1|  <  dQ  (3.71) 

DA>  =  \max{$d,E(<f>d))/min($d,E($d))  -  1|  <  dE.  (3.72) 

After  testing  several  methods  for  satisfying  both  conditions,  we  choose  to  search  for  the  minimum  of  P  =  [DQ  * 
ZA<E»]  (A)  within  the  tolerance  defined  by  dQ  and  dE.  Variation  on  the  definition  of  the  product  P  of  the  two 
conditions  do  not  offer  any  gain. 

With  this  second  method  and  estimation  of  noise  through  the  L-curve  we  find  that  tighter  constraints  are 
applicable.  The  method  is  validated  by  inversion  of  synthetic  data  (produced  with  dipole  model)  and  real  data 
obtained  from  the  USACE-ERDC  test  stand  for  40  mm,  60  mm,  90  mm  and  M42  projectiles,  for  which  two  depths 
(1  shallow,  1  deep)  and  three  inclinations  (0,  45,  90  degrees)  were  surveyed  with  a  26-time-channel  Geonics  EM- 
63.  In  the  case  of  real  data,  the  corner  of  the  L-curve  proved  to  be  an  adequate  solution  in  19%  of  inversions 
with  DQ<12.5%;  for  the  remaining  inversions,  38%  of  cases  were  within  the  bound  dQ=dE=12.5%,  30%  within 
dQ=dE=25%,  12%  within  dQ=2dE=100%  and  the  last  1%  within  dQ=200%  and  dE=75%. 

Application  of  this  method  to  the  40  mm  projectile  is  illustrated  in  Figures  3.27  and  3.28  by  the  large  dots 
previously  defined  as  alternative  regularization.  Extensive  tests  presented  in  the  following  section  confirm  the 
robustness  of  the  method  in  a  wider  setting.  Figure  3.29  illustrates  the  effect  of  different  regularizations  for 
the  inversion  of  noisy  data  acquired  over  a  40  mm  UXO  at  the  Sky  Research  test  plot.  Comparison  between 
observed  data  and  predictions  shows  that  all  models  except  for  the  most  regularized  one  yield  acceptable  fit  with 
the  observed  data.  Differences  between  the  models  for  the  L-curve  and  the  chosen  regularization  are  negligible. 


Conclusion 

We  have  now  established  a  procedure  for  regularized  inversion  of  electromagnetic  data  with  SMC  that  pro¬ 
vides  a  forward  model  that  closely  predicts  surface  observations  and  a  stable  total  magnetic  charge  derived  from 


69 


L-curve 


Corner  L-curve 


Chosen  regularization 


Figure  3.29:  Regularization  of  field  data  over  a  40  mm  UXO.  L-curve,  Observation  and  model  for  different  levels 
of  regularization.  Corner  of  L-curve  indicated  with  a  star,  our  chosen  regularization  with  a  red  dot. 
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the  uniform  charge  density.  The  total  magnetic  charge  derived  with  this  method  is  independent  of  the  number 
of  rings  used  for  discretization,  would  they  be  10,  40  or  60.  In  the  following  applications  40  rings  are  used  to 
allow  for  the  complexity  of  bodies  with  composite  structure  (geometry,  material  with  different  conductivity  and 
magnetic  permeability),  while  regularization  handles  the  redundancy  of  parameters. 

Alternative  avenues  of  regularization  were  also  explored  during  our  close  examination  of  the  surface  mag¬ 
netic  charge  model.  We  replaced,  for  instance,  the  search  for  a  smooth  model  by  one  that  exhibits  the  smallest 
total  magnetic  charge  (replacing  the  elements  of  the  regularization  matrix  Wm  by  the  surface  area  of  each  ring). 
The  resulting  regularization  remained  too  weak  to  temper  out  the  large  oscillations  in  charge  amplitudes.  We  also 
tested  WTO  as  the  second  derivative  operator,  Wm  becoming  either  a  rank-2  deficient  operator,  or  a  full-rank 
matrix  when  boundary  conditions  such  as  null  charge  beyond  the  spheroid  were  applied.  These  formulations  lead 
to  similar  results  to  those  reported  for  the  standard  case.  Black-box  types  of  regularization  procedures  such  as  the 
MatLab  functions  lsqnonneg  (positivity  constraint)  and  linprog  (linear  programming  solutions)  were  found  to  be 
less  stable  than  the  regularizations  discussed  above,  both  for  rings  and  patches.  Another  possible  way  of  control¬ 
ling  the  charge  distribution  is  decomposition  through  a  given  number  of  simple  basis  functions  (sine-cosine).  This 
method  was  tested  on  a  sphere  (icosahedron)  with  two-dimensional  Fourier  series  truncated  to  limited  modes  to 
restrict  the  amount  of  features  in  the  charge  distribution.  Results  not  presented  here  show  that  stability  of  the  total 
magnetic  charge  can  be  achieved  with  careful  choice  of  truncation.  A  systematic  procedure  remains  to  be  formu¬ 
lated.  Because  modeling  an  ordnance  with  a  sphere  simplifies  the  inversion  problem  by  removing  the  orientation 
variables,  further  work  is  warranted  to  explore  that  possible  avenue. 


3.6.4  Stability  of  the  regularization  method 


Adding  controlled  noise 

Ability  of  the  regularized  inversion  to  produce  a  model  with  misfit  || Zm  —  d||  comparable  to  noise  can  be 
tested  in  an  environment  where  noise  is  controlled.  This  can  be  done  by  considering  shallow  measurements  (Hsc) 
with  the  40  mm  projectile  where  signal  is  high  and  noise  remains  low.  Noise  can  then  be  arbitrarily  added  to  the 
instrument  (ejnst,  %)  and  to  the  data  (edata,  %)  according  to 

-  added  instrument  error:  IE  =  R1  *  einst 

-  added  data  error:  DE  =  Hsc  *  diag{R2 )  *  edata 

-  total  added  noise  A E  =  IE  +  DE, 

with  Rl  and  R2  vectors  of  length  M  with  random  entries,  chosen  from  a  normal  distribution  with  mean  zero, 
variance  one  and  standard  deviation  one;  diag(R2)  is  the  matrix  with  diagonal  R2  and  zeros  everywhere  else. 
Results  are  summarized  in  Table  3.3  and  3.4. 


^inst?  ^data 

l|i/wd|| 

HASH 

S($d) 

Q 

^uni 

Corr. 

No  noise 

28 

80 

0 

0.35 

-0.0242 

-0.0260 

1 

500,  10 

100 

84 

95 

1.0 

-0.0248 

-0.0265 

0.99 

2000,  30 

323 

73 

315 

3.8 

-0.0252 

-0.0253 

0.92 

2000,  50 

419 

73 

496 

5.2 

-0.0259 

-0.0273 

0.87 

Table  3.3:  40  mm  object  in  horizontal  position,  25  cm  below  sensor,  1st  time  channel.  <l\i  is  model  misfit; 
1 1  /W,/||  is  the  assigned  error  (from  automated  estimate  of  standard  deviation  of  data  and  instrument  error); 
||AS|  is  noise  added  to  noise-free  data;  A(4yi )  is  the  expected  noise  as  in  W,/,  departure  from  1  means  that 
automated  guess  on  noise  level  was  inaccurate.  The  method  is  deemed  successful  if  <l\i  is  close  to  |  A  /t  |  and  Q 

tO  Quni- 


Both  tables  confirm  that  even  when  data  are  largely  corrupted  by  noise  the  proposed  regularized  inversion 
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^inst?  ^data 

l|i/wd|| 

||AF;|| 

S($  d) 

Q 

Quni 

Corr. 

No  noise 

82 

102 

0 

0.80 

-0.0253 

-0.0254 

1.00 

500,  10 

153 

92 

158 

1.5 

-0.0237 

-0.0239 

0.98 

2000,  30 

442 

99 

403 

4.1 

-0.0242 

-0.0246 

0.89 

2000,  50 

715 

134 

807 

4.84 

-0.0413 

-0.0277 

0.87 

Table  3.4:  40  mm  object,  45°  inclination,  25  cm  below  sensor,  1st  time  channel. 


method  maintains  (1)  stable  total  magnetic  charge,  (2)  misfit  comparable  with  added  noise  ||AS|  and  (3)  large 
correlation  between  data  and  prediction.  The  third  columns  in  both  tables  suggest  that  the  method  for  assigning 
expected  errors  in  the  weighting  matrix  Wd  is  not  optimal,  expected  errors  are  generally  over-estimated  when 
noise  is  low  and  under-estimated  when  noise  is  high.  This  shortcoming  does  not,  however,  affect  the  performance 
of  the  proposed  regularization  method,  which  relies  instead  on  the  corner  of  the  L-curve  to  identify  eligible 
models  that  strike  the  right  balance  between  model  complexity  and  misfit.  This  property  is  confirmed  by  the 
general  agreement  between  added  noise  and  misfit  with  the  chosen  model.  It  is  also  confirmed  that  amongst  those 
eligible  models  it  is  possible  to  select  some  with  total  magnetic  charge  near  Qunj.  Results  hold  for  large  amounts 
of  noise:  relative  error  for  Q  reaches  62%  for  added  noise  ejata  =  50%  in  Table  3.5. 


Noise  and  regularization  for  different  types  of  ordnance 

Adaptability  of  the  method  to  different  size  and  shape  of  ordnance  is  tested  on  three  additional  standard 
items:  the  60  mm,  90  mm  and  M42  ordnance.  Several  orientations  and  depths  are  tried  (Table  3.5)  with  the  1st 
and  12th  time  channels  of  data  collected  with  a  Geonics  EM-63  sensor  at  the  USACE-ERDC  Vicksburg,  MS  test 
stand.  These  items  are  chosen  because  the  40  mm,  60  mm  and  90  mm  have  distinct  sizes,  while  the  40  mm  and 
M42  have  similar  size  but  different  physical  properties  (40  mm  is  made  out  of  aluminum  and  copper,  M42  of 
steel).  The  1st  time  channel  is  taken  for  its  high  SNR  and  the  later  12th  time  channel  to  illustrate  the  effect  of 
high  noise  levels.  Figure  3.30  and  3.31  show  results  of  inversions  for  data  processed  with  their  original  noise. 
Table  3.6  gathers  significant  statistics  for  the  total  magnetic  charge  derived  from  the  uniform  distribution. 


Object  /  Inclination 
Depth 

Horiz. 

Shallow 

45° 

Shallow 

Vert. 

Shallow 

Horiz.,  45°,  Vert. 
Deep 

40  mm 

25.4  cm 

21  cm 

24  cm 

60  cm 

60  mm 

60  cm 

60  cm 

60  cm 

100  cm 

90  mm 

60  cm 

60  cm 

60  cm 

100  cm 

M42 

25.8  cm 

21.8  cm 

25  cm 

60  cm 

Table  3.5:  Depth  and  inclination  of  measurements  over  test  stand. 

Figure  3.30  shows  results  for  the  1st  time  channel.  In  all  panels  the  total  magnetic  charge  exhibits  large 
variations  during  regularization,  except  when  the  item  is  not  horizontal,  as  previously  observed  in  Figure  3.28. 
Choice  of  the  corner  of  L-curve  (star  markers)  yields  unstable  total  magnetic  charge,  especially  for  60  mm  and 
M42  items.  Conversely,  Table  3.6  show  that  Quni  obtained  with  the  uniform  charge  model  keeps  a  stable  value 
through  changes  of  position,  its  standard  deviation  remaining  at  least  five  times  smaller  than  the  amplitude  of 
Quni.  Use  of  regularization  for  Q  near  Quni  (dot  markers  in  Figure  3.30)  provides  both  stable  total  magnetic 
charge  and  adequate  fit  to  data. 

Figure  3.31  shows  the  same  four  items  at  the  12th  time  channel,  where  signal  gets  weaker  and  noise  relatively 
increases.  Examination  of  SNR  reveals  that  at  a  shallow  depth  clear  distinction  of  the  target  is  possible  with  SNR 
larger  than  12.  At  deeper  positions  inversion  of  data  remains  possible  even  though  SNR  lowers  to  2-8  (minimum 
for  60  mm  item  horizontal  at  1  m  below  sensor  and  M42  horizontal  at  60  cm  below  sensor,  both  shown  with 
magenta  line  with  plus  signs).  Study  of  late  time  channels  and  noisy  environments  suggests  increasing  separation 
between  “classic”  (L-curve)  and  “special”  (L-curve+Quni)  regularization  methods  when  SNR  decreases,  here 
visible  for  all  four  items  in  Figure  3.31.  With  different  choices  of  the  regularization  parameter  Q  can  vary  tenfold 
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Figure  3.30:  Variation  of  Total  magnetic  charge  Qsmc  during  regularized  inversion  of  first  time  channel  for  60  mm, 
90  mm  and  M42  projectiles.  For  all  panels,  each  curve  shows  Qsmc  for  a  different  position/orientation  of  target 
relative  to  test  stand.  Horizontal  dashed  lines  show  Q uni  for  each  case  and  define  the  region  of  stability  of  the  total 
charge.  Large  dots  indicate  Qsmc  for  the  proposed  regularization  choice  whereas  stars  highlight  failure  to  obtain 
stable  charge  with  corner  of  L-curve.  Positions/orientations  are  given  in  Table  3.5.  Black,  red  and  green  curves 
correspond  to  shallow  depths,  blue,  magenta  and  cyan  to  deeper  positions;  black  solid  line  and  blue  line-dot  at  45 
degrees,  red  dashed  and  cyan  dash-dot  vertical,  green  crosses  and  magenta  pluses  horizontal,  as  in  Figure  3.27. 


UXO/Quni 

Mean  Quni(T1) 

&Quni(Tl) 

Mean  Quni(T12) 

crQuni(Tl2) 

40  mm 

-0.023 

0.0024 

-0.0076 

0.0007 

60  mm 

-0.326 

0.062 

-0.073 

0.022 

90  mm 

-0.915 

0.114 

-0.206 

0.051 

M42 

-0.041 

0.007 

-0.003 

0.0009 

Table  3.6:  Mean  and  standard  deviation  of  the  total  magnetic  charge  for  uniform  distribution  of  normalized 
surface  charge  density,  at  1st  (Xj)  and  12th  (7 ’12)  time  channels. 


reg 


reg 


reg 


reg 


Figure  3.31:  Total  magnetic  charge  Qsmc  during  regularization  for  40  mm,  60  mm,  90  mm  and  M42.  Each  panel 
shows  all  six  positions  given  in  Table  3  using  the  12th  time  channel  with  the  same  plotting  conventions  as  previous 
figure. 
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and  even  change  sign,  as  observed  for  the  60  mm  and  M42,  while  Qunj  remains  stable  for  all  items  (Table  3.6). 


3.6.5  Detailed  analysis  of  four  standard  ordnance 


Total  magnetic  charge  and  time 


Figure  3.32:  Total  magnetic  charge  as  a  function  of  time  for  40  mm,  60  mm,  90  mm  and  M42  projectiles,  using 
EM-63  data  from  test  stand  at  positions  summarized  in  Table  3.  Time  series  are  shorter  for  smaller  objects  because 
SNR  is  too  low  at  later  time  channels. 

As  a  last  test  on  method  stability,  regularized  inversion  and  detailed  analysis  are  applied  to  data  sets  for 
40  mm,  60  mm,  90  mm  and  M42  standard  ordnance  at  multiple  depths  and  orientations.  The  Geonics  EM-63 
sensor  acquired  data  on  the  USACE-ERDC  Vicksburg  test  stand  at  26  time  channels  (0.18  ms  to  25.14  ms),  two 
depths  (one  shallow,  one  deep)  and  three  inclinations  (0,  45,  90  degrees),  as  in  Table  3.5.  Figure  3.32  shows  the 
total  magnetic  charge  as  a  function  of  time  for  all  measured  configurations  obtained  with  the  proposed  regularized 
inversion  algorithm  for  NSMC.  Several  conclusions  can  be  drawn: 


•  All  lines  cluster  for  each  object;  the  recovered  total  magnetic  charge  is  therefore  a  stable  feature  and  the 
inversion  is  robust; 

•  Each  object  has  a  distinct  time  evolving  total  magnetic  charge; 

•  The  magnitude  of  the  total  magnetic  charge  scales  with  the  volume  of  the  object; 

•  Objects  with  different  physical  properties  have  different  time  decay  for  the  total  magnetic  charge,  as  illus¬ 
trated  by  the  M42  and  40  mm  items  that  have  similar  size  but  different  material. 
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These  results  suggest  clear  and  stable  separation  with  NSMC  and  the  total  magnetic  charge  for  these  four  types 
of  ordnance,  thus  opening  the  possibility  of  applying  automated  discrimination  procedures.  As  a  side  note,  there 
is  a  slight  increase  of  the  total  magnetic  charge  for  the  10th  and  11th  time  channels.  This  effect  is  only  due  to  a 
pervasive  instrument  bias,  not  to  any  physical  property  or  modeling  issue,  and  would  not  appear  should  the  sensor 
be  perfectly  calibrated. 


Error  analysis 


Inclination  of  UXO  (°) 

Figure  3.33:  Correlation  coefficient  between  observed  data  and  model  prediction  for  all  time  channels,  2  depths 
and  3  inclinations  (0  for  horizontal). 

Having  established  the  stability  of  the  total  magnetic  charge  for  the  40  mm,  60  mm,  90  mm  and  M42  at 
various  depths  and  orientations,  we  examine  the  quality  of  prediction  of  the  forward  models  employed  above. 
Special  consideration  is  given  to  correlation,  misfit  and  ability  to  perform  at  low  SNR.  Figure  3.33  shows  the 
correlation  coefficient  between  the  observation  and  model  prediction  for  all  time  channels  as  a  function  of  the 
inclination  of  the  item,  with  the  different  depths  superposed.  Correlation  close  to  unity  in  most  cases  proves  that 
model  predictions  are  reliable.  The  few  occurrences  of  low  correlation  correspond  to  late  time  channels  where 
data  quality  is  poor.  Best  performance  is  achieved  for  targets  in  vertical  position,  in  which  case  the  signal  recorded 
at  the  surface  takes  the  shape  of  a  simple  radial  decay  away  from  the  center  of  target  because  UXO  are  bodies  of 
revolution.  Success  when  there  is  central  symmetry  in  the  data  is  not  surprising  because  our  model  does  assume 
central  symmetry  of  charges  along  the  main  axis  of  the  modeling  spheroid  (i.e.,  rings  with  uniform  charge)  here 
in  vertical  position,  and  thus  naturally  fits  the  data. 

Figure  3.34  measures  the  misfit  of  models  obtained  through  the  proposed  regularized  inversion  method  with 
the  expected  noise  estimated  from  the  corner  of  the  L-curve.  In  97%  of  tested  cases  misfit  lies  within  a  25% 
relative  difference  with  estimated  noise,  thus  proving  that  the  method  generally  converges. 

Figure  3.35  provides  the  conditions  of  signal  to  noise  ratio  throughout  the  sets  of  measurements.  SNR  varies 
by  several  orders  of  magnitude  for  different  depths  and  times.  Inversions  are  not  carried  out  when  SNR  is  lower 
than  two  or  when  there  are  less  than  ten  data  with  values  above  the  estimated  standard  deviation  of  noise;  SNR  is 
then  set  to  0.  SNR  is  particularly  low  at  late  time  channels,  or  at  all  times  for  small  items  placed  at  greater  depth, 
e.g.  40  mm,  60  mm  and  M42.  Regularized  inversion  of  NSMC  is  able  to  perform  at  low  SNR  and  has  potential 
for  use  even  with  highly  contaminated  data. 
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5  10  15  20  25 


Time  channel 

Figure  3.34:  Percentage  of  solutions  for  which  ^{regu)  /  §  ^(corner)  —  1|  smaller  than  misfit  ratio  (79$)  as  a 
function  of  time  channel.  The  lowest  row  shows  that  in  55%  of  cases  the  misfit  ratio  is  below  5%,  which  means 
that  the  misfit  for  the  chosen  parameter  is  within  5%  of  that  of  the  corner  of  the  L-curve  (assumed  best  guess  for 
noise  estimation).  The  next  lines  show  that  73%  of  cases  are  below  10%  misfit  ratio,  97%  below  25%,  99.5% 
below  50%  and  100%  below  100%. 
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Figure  3.35:  Square  root  of  signal  to  noise  ratio.  For  each  item,  the  first  three  columns  correspond  to  shallow 
positions,  the  next  three  to  deeper  positions,  where  SNR  is  lower.  SNR=0  at  later  time  channels  if  too  few  data 
for  inversion  (22%  of  cases),  SNR>2  otherwise.  Inversion  was  successful  even  for  SNR  as  low  as  2. 
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3.7  Application  of  the  regularized  NSMC 

3.7.1  Time  domain  analysis  of  standard  items 


Fifteen  standard  ordnance  were  measured  at  the  USACE-ERDC  Vicksburg  test  stand  using  a  Geonics  EM-63  time 
domain  sensor  to  test  and  develop  the  capabilities  of  current  EMI  forward  models.  These  ordnance  span  a  wide 
range  of  sizes,  from  155  mm  to  20  mm  projectiles,  and  different  material.  Data  were  taken  at  two  depths  and  three 
orientations,  as  in  the  previous  section.  Figure  3.36  shows  total  magnetic  charge  for  all  these  items.  To  avoid 
saturating  the  figure  only  the  median  among  all  positions  at  each  time  channel  is  displayed.  Not  shown  here,  the 
distribution  of  total  magnetic  charge  with  time  when  depth  and  orientation  has  similar  standard  deviation  as  the 
four  items  presented  in  the  previous  section.  This  illustrates  that  the  method  of  deriving  the  total  magnetic  charge 
is  stable  for  all  types  of  ordnance  at  hand. 

Several  observations  can  be  made.  The  amplitude  of  total  magnetic  charge  decays  with  time  as  the  amplitude 
of  the  scattered  field  does.  At  early  time,  the  magnitude  of  the  total  magnetic  charge  scales  with  the  size  of  the 
ordnance.  At  later  times,  the  charge  signal  reflects  other  properties  of  ordnance  (shape,  material).  The  latest  time 
channels  are  not  inverted  for  small  items  because  the  signal  is  too  weak.  The  forward  model  provides  adequate 
fit  to  data  and  high  degrees  of  correlation  between  observations  and  predictions  are  achieved.  All  data  were 
successfully  inverted  on  the  same  size  of  spheroid,  a  priori  information  on  the  size  of  the  object  to  model  is 
therefore  not  required. 

These  results  indicate  that  the  time-evolving  total  magnetic  charge  acts  as  a  unique  signature  for  each  type 
of  ordnance.  For  instance  the  37  mm,  MK118  and  BLU26,  which  have  similar  diameters  but  different  material, 
length  and  density,  have  similar  amplitudes  at  early  time  channels  whereas  their  relative  charge  varies  with  several 
orders  of  magnitude  beyond  the  15th  time  channel.  A  preliminary  attempt  to  quantify  the  degree  of  separation 
between  each  class  of  ordnance  is  presented  in  Figure  3.37  in  a  simple  canonical  analysis  as  in  Beran  (2005), 
using  the  first  18  time  channels.  In  the  left  panel  most  of  the  variability  recorded  amongst  all  test  stand  data  for 
the  15  ordnance  can  be  summarized  with  two  eigenvalues  and  feature  vectors  xc\  and  c2.  The  resulting  two- 
dimensional  space  xc\  —  xc‘2  (right  panel)  of  Figure  3.37  has  total  charges  corresponding  to  different  positions 
for  each  ordnance  cluster  and  separate  from  those  of  other  ordnance.  The  NSMC  appears  to  be  suitable  for 
discrimination  algorithms.  Further  work  is  warranted  to  extract  the  most  relevant  information  out  of  the  total 
charge  distributions  and  to  deal  with  complex  situations,  such  as  discrimination  with  incomplete  time  series  of 
total  charge. 

As  a  further  discrimination  test,  a  series  of  cylinders  were  measured  over  the  USACE-ERDC  test  stand  in  the 
same  conditions  as  the  standard  ordnance  to  study  their  EMI  response  and  test  models.  These  six  cylinders  were 
chosen  to  be  solid  or  hollow,  short  or  long,  made  of  steel  or  aluminum.  Their  total  magnetic  charge  is  presented 
in  Figure  3.38,  along  with  three  large  fragments  of  105  mm  HEAT  rounds.  Comparison  with  Figure  3.36  shows 
that  the  charge  of  these  cylinders  has  a  smaller  time  decay  rate,  which  distinguishes  the  cylinders  from  their  UXO 
counterparts.  The  fragments  of  105  mm  also  produce  distinct  total  charge. 


3.7.2  Frequency  domain  analysis  of  standard  items 

Data  were  also  collected  with  a  Geophex  GEM-3  sensor  at  the  USACE-ERDC  Vicksburg  test  stand  to  test  the 
models’  ability  to  perform  with  frequency-domain  electromagnetic  data.  The  sensor  provides  the  real  and  imag¬ 
inary  part  of  the  signal  sampled  at  10  frequencies  from  90  to  41010  Hz.  Data  were  processed  and  inverted  by 
taking  the  modulus  of  data  for  every  time  channel  so  that  amplitude  reached  maximum  value  directly  above  target, 
easing  comparison  with  uniform  charge  distribution  model  so  that  the  exact  same  regularized  inversion  procedure 
as  for  time  domain  could  be  applied.  Information  on  phase  was  not  used.  Results  are  presented  in  Figure  3.39. 

The  regularized  inversion  of  NSMC  performs  as  well  in  the  frequency  domain  as  in  the  time  domain.  Pre¬ 
dicted  total  magnetic  charge  is  consistent  for  different  positions  and  orientations  of  the  item.  The  data  misfit  is 
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Figure  3.36:  Total  magnetic  charge  of  15  standard  ordnance  as  a  function  of  time.  Data  acquired  over  test  stand 
with  EM-63  sensor  at  26  channels.  Only  median  is  shown  for  clarity,  standard  deviations  similar  to  those  of 
Figure  3.32.  Note:  there  are  2  types  of  81  mm  items,  the  ATC  and  Montana. 


78 


Figure  3.37:  Canonical  analysis  of  the  15  standard  UXOs  shown  in  Figure  3.36.  Left  panel:  eigenvalues.  Right 
panel:  representation  of  features  from  total  magnetic  charges  taken  at  different  depths  and  orientations  in  plane  of 
the  two  first  eigenvectors. 


close  to  expected  noise  and  correlation  coefficients  are  close  to  unity,  but  those  details  are  not  presented  here. 
The  method  is  therefore  robust  and  the  magnitude  of  the  total  charge  scales  with  the  size  of  objects.  Canonical 
analysis  presented  in  Figure  3.40  confirms  stability  of  the  total  magnetic  charge  obtained  by  regularized  inversion, 
and  shows  that  each  class  of  UXO  has  a  distinct  signature  and  belongs  to  a  different  part  of  feature  space.  These 
results  show  similar  performance  for  regularized  inversion  of  the  surface  charge  model  can  indeed  be  expected 
for  time  domain  and  frequency  domain  electromagnetic  induction-based  identification  of  ordnance. 


3.7.3  Application  to  field  data 


The  following  presents  a  series  of  preliminary  tests  to  assess  the  type  of  difficulties  the  regularized  inversion  of 
NSMC  would  face  when  confronted  with  field  realities. 


Discrimination  by  total  magnetic  charge 

Most  of  our  development  work  on  the  regularized  inversion  of  NSMC  presented  so  far  is  based  on  test 
stand  data,  which  offers  a  mostly  noise-free  environment.  As  a  first  step  toward  dealing  with  real  site  conditions, 
inversion  was  performed  on  data  collected  over  a  test  plot  at  the  Sky  Research  center  in  Ashland,  Oregon,  where 
a  series  of  standard  UXO  were  seeded.  Data  were  acquired  with  an  EM-63  sensor  in  cued  interrogation  mode, 
that  is  a  stable  platform  with  fix  stations,  as  opposed  to  the  regular  dynamic  acquisition  mode  of  real  survey.  Data 
used  here  were  taken  40  cm  above  the  surface,  slightly  higher  than  the  position  of  the  sensor  on  its  cart,  which 
partially  reduces  the  electromagnetic  response  of  the  magnetic  soil.  The  following  analysis  includes  a  uniform 
background  as  additional  parameter  in  the  inversion  to  solve  for  the  additional  EM  effect  of  soil. 

Inversion  results  for  three  unknown  targets,  for  which  the  position  was  estimated  by  other  means,  are  pre¬ 
sented  in  Figure  3.41.  Thick  lines  show  their  total  magnetic  charge  (derived  from  the  field  data),  thin  lines  recall 
the  total  magnetic  charge  obtained  from  the  test  stand,  as  in  Figure  3.32.  Recovered  total  charges  clearly  overlap 
with  their  test  stand  counterparts;  all  three  items  are  identified  despite  significant  noise  that  corrupts  late  time 
channels.  Results  for  the  40  mm  and  60  mm  ordnance  are  unambiguous  for  all  time  channels,  whereas  the  M42 
drifts  toward  a  40  mm-type  response  at  later  times.  Work  in  progress  is  trying  to  address  that  issue. 
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Figure  3.38:  Total  magnetic  charge  of  six  types  of  cylinders  as  a  function  of  time  using  EM-63  data  taken  over 
a  test  stand.  SSL:  solid  steel  long,  SSS:  solid  steel  short,  HSL:  hollow  steel  long,  HSS:  hollow  steel  short,  SAL: 
solid  aluminum  long,  SAS:  solid  aluminum  short. 
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Figure  3.39:  Total  magnetic  charge  of  13  standard  UXO  items  as  a  function  of  time,  for  data  acquired  with  an 
GEM-3  sensor  over  Vicksburg  test  stand.  Only  median  is  shown  for  clarity,  variance  is  illustrated  in  canonical 
analysis  of  Figure  3.40. 
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Figure  3.40:  Canonical  analysis  on  GEM-3  data  for  13  classes  of  items.  Left  panel:  total  magnetic  charge  obtained 
for  all  positions  measured  on  test  stand.  Central  panel:  eigenvalues.  Right  panel:  representation  of  features  from 
total  magnetic  charges  taken  at  different  depths  and  orientations  in  plane  of  the  two  first  eigenvectors. 


Figure  3.41:  Field  test  of  NSMC.  Total  magnetic  charge  for  three  buried  items,  shown  in  thick  lines,  overlaid  on 
top  of  the  total  magnetic  charge  derived  for  the  same  ordnance  on  test  stand. 
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As  a  side  note,  inversion  of  field  data  proved  to  be  beneficial  for  the  development  of  our  inversion  procedure 
because  it  forced  a  clear  definition  of  the  uniform  distribution.  When  regularization  was  introduced  (Section  3.6), 
the  uniform  charge  distribution  was  taken  as  the  median  of  the  distribution  of  fT|c /HjZij,  the  ratio  of  the  measured 
scattered  field  and  the  sum  of  rows  in  the  modeling  matrix.  The  distribution  and  median  change  when  noise  and 
background  value  increase.  The  most  robust  mean  of  defining  the  uniform  charge  distribution  is  that  of  the  upper 
end  member  of  regularization.  This  definition  of  the  uniform  charge  proves  to  be  far  more  robust  because  the 
regularization  process  removes  the  effect  of  noise.  Further  examination  of  field  data  shall  confirm  the  success  of 
the  method. 


Effect  of  position  on  total  magnetic  charge 

Defining  the  center  of  a  target  for  EMI  modeling  is  not  straightforward  because  ordnance  have  complex 
shapes  and  include  different  material.  In  SMC  modeling,  ordnance  are  represented  by  a  spheroid,  the  center  of 
which  does  not  necessarily  coincide  with  that  chosen  for  test  stand  measurements.  This  leads  to  positional  un¬ 
certainties  for  inversion  of  test  stand  data.  Addressing  that  issue  and  others  regarding  definitions  of  orientation, 
we  conducted  a  simple  sensitivity  analysis  and  found  that  the  total  magnetic  charge  obtained  by  NSMC  showed 
negligible  sensitivity  to  horizontal  positioning  error  within  15  cm  and  to  any  amount  of  error  in  azimuthal  direc¬ 
tion.  The  misfit,  however,  greatly  increased  when  such  errors  occurred  and  therefore  provided  a  good  tool  for 
estimating  the  correct  parameters. 


Time 


Figure  3.42:  Effect  of  depth  on  the  total  magnetic  charge.  Upper  panel:  charge  vs.  time;  lower  panel:  correlation 
between  observation  and  prediction. 

Depth  is  a  more  complex  issue  that  requires  particular  attention.  Figure  3.42  shows  the  effect  of  modeling 
the  wrong  depth  on  the  total  magnetic  charge  of  a  40  mm  UXO.  The  upper  panel  shows  the  charge,  the  lower 
panel  the  correlation  between  the  data  and  forward  model  predictions.  The  total  charge  is  found  to  vary  in  a 
predictable  manner  by  a  factor  two  for  a  10  cm  increment  in  vertical  direction,  expectedly  increasing  upward  and 
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Figure  3.43:  Percent  change  of  the  total  magnetic  charge  when  the  model  depth  is  offset  from  the  true  depth. 

decreasing  downward  with  respect  to  the  correct  location.  Moreover,  the  charge  presents  exactly  the  same  time 
decay  at  all  modeled  depths.  Given  the  large  separation  between  the  total  charge  and  the  distinct  decay  rates  of 
each  UXO,  this  result  suggests  that  10-20  cm  errors  in  depth  would  be  no  impediment  for  identifying  the  buried 
item.  From  there,  recognizing  the  buried  item  can  add  a  constraint  on  the  inversion  to  recover  the  correct  burial 
depth.  The  lower  panel  shows  that  the  correlation  coefficient,  which  reflects  the  data  misfit,  is  only  moderately 
but  consistently  affected  by  the  depth  error. 

Figure  3.43  shows  a  first  attempt  to  generalize  that  study.  EM-63  test  stand  data  are  inverted  adding  an  offset 
to  the  true  depth  (vertical  axis),  different  settings  of  depth  and  orientation  for  each  targets  are  indicated  on  the 
horizontal  axis,  only  the  first  time  channel  is  presented.  The  total  charge  varies  by  50-100%  within  10  cm  of  the 
true  position,  which  confirms  results  for  the  40  mm  ordnance.  Detailed  examination  through  formal  sensitivity 
analysis  is  warranted  to  reinforce  that  preliminary  result.  Confirmation  of  that  effect  would  make  the  NSMC 
method  a  potent  tool  for  UXO  identification  and  discrimination. 


3.8  Conclusion 


In  this  chapter  we  conducted  a  detailed  examination  of  the  normalized  surface  magnetic  charge  model  to  apply  it 
to  electromagnetic  data.  With  the  goal  of  using  the  total  magnetic  charge,  the  integral  of  surface  magnetic  charges, 
as  a  criterion  to  identify  and  discriminate  buried  ordnance  against  clutter,  we  particularly  focused  on  the  stability 
of  TMC  for  multiple  orientations,  depths,  times  and  sensor  types.  Because  the  inverse  problem  is  ill-posed,  direct 
inversion  is  unstable  and  requires  regularization.  We  looked  for  a  physical  basis  to  carry  on  that  regularization  and 
investigated  the  properties  of  the  normalized  surface  charge  distribution  at  the  surface  of  a  spheroid  enclosing  an 
ordnance.  We  found  that  the  distribution  varied  under  different  directions  and  positions  of  excitation.  When  data 
are  collected  in  real  life,  targets  are  illuminated  under  several  directions,  therefore  some  sort  of  averaging  of  the 
charge  distribution  occurs.  This  theoretical  consideration,  coupled  with  extensive  empirical  analysis  of  synthetic 
and  test  stand  data  taken  at  multiple  depths  and  angles,  leads  to  the  conclusion  that  the  most  stable  feature  of  the 
surface  charge  distribution  and  its  total  charge  is  the  uniform  charge  distribution  obtained  by  averaging.  This  result 
forms  the  basis  of  a  regularization  procedure  that  identifies  a  forward  model  that  simultaneously  (1)  produces  a 
high  correlation  coefficient  between  observation  and  prediction,  (2)  balances  the  misfit  (not  fitting  noise)  and  the 
complexity  of  the  charge  distribution,  and  (3)  produces  a  total  magnetic  charge  close  to  that  of  a  uniform  charge 
distribution. 

The  method  was  tested  on  data  collected  with  the  Geonics  EM-63  and  Geophex  GEM-3  sensors  over  15 
standard  ordnance,  6  types  of  cylinders  and  clutter  at  multiple  depths  and  orientations  with  the  help  of  USACE 
at  the  ERDC  test  stand  in  Vicksburg,  MS.  We  proved  that  the  total  magnetic  charge  derived  from  this  regularized 
inversion  protocol  was  a  robust  and  stable  feature  for  each  type  of  target  with  both  time-domain  and  frequency- 
domain  types  of  data.  We  also  found  that  the  amplitude  and  time  decay  of  TMC  differs  for  each  class  of  ordnance 
and  cylinders  and  reflects  the  size  and  electromagnetic  properties  of  their  material.  A  simple  canonical  analysis 
confirmed  the  clear  separation  between  the  TMC  features  of  all  ordnance. 

Preliminary  tests  were  also  done  to  assess  the  robustness  of  the  method  with  noisy  field  data  and  the  sen- 
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sitivity  of  target  positioning  error.  Inversion  of  field  data  taken  at  the  Sky  Research  test  plot  was  successful  at 
identifying  buried  items  despite  the  high  noise  generated  by  the  magnetic  soil.  Sensitivity  analysis  when  the  target 
position  and  orientation  were  altered  showed  that  large  errors  in  horizontal  position  and  azimuthal  direction  do  not 
affect  determination  of  TMC.  In  contrast,  the  amplitude  of  TMC  scales  with  depth  error  while  time  decay  remains 
unaffected,  thus  allowing  possible  identification  of  the  target  despite  depth  error.  Further  work  to  systematically 
investigate  these  effects  and  others  is  warranted  and  detailed  in  the  following  Conclusion  chapter. 
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Chapter  4 


Discussion  and  Conclusions 

4.1  Objectives 


Initially,  the  objective  of  this  research  project  was  to  apply  the  Method  of  Auxiliary  Sources  as  a  modeling  tool 
for  inversion  of  UXO  sensor  data.  Through  comparisons  between  measured  data,  analytical  solutions,  and  data 
predicted  by  MAS,  we  found  that  MAS  was  able  to  accurately  model  the  EMI  response  of  compact  metallic 
objects  such  as  UXO  over  a  large  frequency  range  extending  to  the  static  case  (Beran,  2005,  Chapter  2).  The 
ability  to  model  over  such  a  wide  range  made  the  MAS  an  attractive  modeling  technique  for  joint  inversion  of 
multiple  data  sets.  Practical  considerations  such  as  the  computational  time  required  by  MAS  to  calculate  sensor 
data  prevented  us  from  applying  MAS  for  inversion  purposes.  However,  its  quality  as  a  physical  model  lead 
us  to  revise  the  objective  of  this  study  to  evaluate  the  potential  of  two  modeling  techniques  derived  from  the 
MAS  framework:  the  Standardized  Excitation  Approach  and  the  Surface  Magnetic  Charge  model.  We  tested  the 
methods  for  their  performance  as  forward  modeling  technique,  assessed  as  their  ability  to  predict  the  data,  provide 
stable  feature  vectors  that  are  unique  to  the  target  and  be  fast  in  the  perspective  of  field  application.  Although  both 
SEA  and  SMC  are  based  on  magnetic  charge  distributions,  they  represent  two  fundamentally  different  approaches 
to  UXO  discrimination.  We  envision  the  SEA  as  part  of  a  library  or  hypothesis  testing  technique,  whereas  SMC 
would  be  used  in  a  parameter  estimation/statistical  classification  scheme. 


4.2  The  Standardized  Excitation  Approach 


In  the  SEA,  the  fundamental  step  is  to  build  up  an  RSS  that  represent  target  attributes  (size,  shape,  and  EM  param¬ 
eters).  We  employ  the  MAS  to  generate  the  RSS  for  a  candidate  UXO  whose  geometry  and  physical  quantities 
are  the  inputs  for  the  process.  The  geometrical  parameters  are  obtained  by  either  digitizing  or  approximating 
equivalent  cylindrical  sections  for  a  composite  object.  The  conductivity  and  magnetic  permeability  of  each  sec¬ 
tion  are  estimated  using  measured  data.  The  MAS  code  is  then  used  to  determine  the  source  distribution  for 
different  spheroidal  modes.  There  is  significant  computational  effort  required  to  generate  the  RSS  and  therefore 
it  cannot  be  used  in  a  parameter  estimation  sense,  i.e.,  it  is  not  reasonable  to  invert  sensor  data  for  an  RSS  (which 
would  subsequently  be  part  of  a  feature  vector  to  be  input  to  a  classification  algorithm).  However,  once  the  RSS 
is  generated,  the  computational  times  are  relatively  quick.  The  SEA/RSS  approach  lends  itself  to  a  library  or  hy¬ 
pothesis  testing  technique  for  discrimination.  Discrimination  is  thus  achieved  by  determining  which  target  within 
the  library  has  the  greatest  likelihood  in  producing  the  anomaly. 

For  this  study  we  developed  an  RSS  library  for  9  UXO.  Testing  was  completed  to  determine  the  number  of 
spheroidal  modes  sufficient  for  modeling  transmitter  fields  generated  by  loop  sensors.  The  ability  to  model  both 
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frequency  and  time  domain  data  was  confirmed  by  comparing  modeling  results  to  Geophex  GEM3  and  Geonics 
EM63  data  collected  on  a  test  stand.  A  simple  library  based  discrimination  technique  was  tested  on  GEM3 
test  stand  data.  An  inversion  algorithm  for  determining  optimal  location  and  orientation  for  a  given  UXO  is  in 
development. 

Some  SEA  issues  were  identified  during  our  investigation: 

•  Improving  the  method  of  estimating  conductivity  and  permeability  for  a  composite  UXO 

When  developing  the  RSS,  the  conductivity  and  permeability  were  calculated  from  a  single  sounding  of 
GEM3  data  collected  over  a  horizontal  target.  A  trial  and  error  procedure  was  used  to  determine  perme¬ 
ability  and  conductivity.  Data  collected  on  a  line  along  the  length  of  the  target  were  used  for  evaluating  the 
suitability  of  the  permeability  and  conductivity  estimates.  Firstly,  the  material  properties  for  an  heteroge¬ 
neous  target  should  not  be  determined  from  only  a  single  sounding  when  additional  spatial  data  is  available, 
or  could  be  collected.  Secondly,  this  process  should  be  reformulated  as  an  inverse  problem  to  determine  the 
material  properties  in  an  optimal  way. 

•  Further  testing  is  required  to  determine  the  optimal  number  of  charge  rings 

We  chose  12  rings  for  the  RSS  based  on  limited  measurements.  Now  that  we  have  access  to  data  acquired 
on  the  USACE-ERDC  test  stand,  we  can  determine  the  number  of  rings  required  by  the  method. 

•  Optimizing  code  for  reducing  forward  modeling  times 

The  non-optimized,  research  code  used  in  this  report  required  close  to  7  seconds  to  model  a  RSS  for  a 
single  sounding  of  13  frequencies.  These  computations  were  carried  out  on  a  Pentium  4  3.2  GHz  processor. 
Preliminary  code  improvements  have  reduced  the  time  to  just  under  three  seconds.  Further  inefficiencies  in 
the  code  are  yet  to  be  fixed.  Eventually,  The  real  high-speed  SEA  computation  might  be  implemented  via 
parallelization. 

•  Developing  algorithms  for  inverting  location  and  orientation 

The  method  tested  in  this  report  determined  location  and  orientation  by  simply  trying  numerous  depths  and 
orientations.  Although  code  was  developed  for  inverting  location  and  orientation  for  a  given  RSS,  it  was 
not  mature  enough  by  the  end  of  the  study  to  evaluate  its  accuracy  and  feasibility  for  real  data. 


4.3  The  Surface  Magnetic  Charge  Model 


The  SMC  approach  is  used  to  predict  the  EMI  response  of  a  metallic  object  by  assuming  that  the  scattered  field 
measured  by  sensors  at  the  surface  originates  from  a  magnetic  charge  distribution  on  a  fictitious  spheroid  that 
encloses  the  target.  Once  the  surface  is  determined,  the  modeled  data  is  a  linear  function  of  the  charge  distribu¬ 
tion  on  the  surface.  The  forward  modeling  is  very  fast  and,  therefore,  sensor  data  can  be  inverted  for  the  charge 
distribution.  The  sum  of  the  charge  distribution,  i.e.,  the  total  surface  magnetic  charge,  can  be  used  for  discrimi¬ 
nation.  Difficulties  arise,  however,  because  determining  the  surface  magnetic  charge  distribution  is  an  ill-posed, 
under-determined  problem.  The  SMC  inversion  problem  requires  significant  regularization  in  order  to  produce 
stable  results  for  the  total  surface  magnetic  charge.  Great  effort  was  invested  in  the  reported  study  to  gain  better 
understanding  of  the  magnetic  charge  distribution  and  accordingly  device  a  stable  regularization  procedure. 

Some  SMC  specific  issues  were  identified  during  our  investigation: 


•  The  Normalized  Magnetic  Charge  distribution  is  dependent  on  the  angle  of  excitation 

We  showed  that  the  normalized  charge  distribution  changes  as  the  angle  of  the  primary  field  changes.  Con¬ 
sequently,  the  total  normalized  magnetic  charge,  suggested  by  Shubitidze  et  al.  (2005a)  as  a  discriminant, 
will  also  vary.  In  order  for  the  total  normalized  magnetic  charge  to  be  a  more  stable  discriminant,  data 
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should  be  taken  such  that  there  are  numerous  excitation  angles.  Data  collected  in  this  manner  would  pro¬ 
duce  a  total  normalized  magnetic  charge  that  represents  an  ’’averaging”  of  the  normalized  charges  over  all 
excitation  directions.  This  type  of  data  could  be  collected  by  a  single  transmitter  loop  in  low  noise  con¬ 
ditions  collecting  multiple  soundings  (such  as  a  test  stand),  or  by  ’’steering”  the  transmit  field  by  using 
combinations  of  transmitter  loops,  or  by  changing  the  orientation  of  a  single  loop. 

•  Regularized  inversion  for  a  stable  Total  Magnetic  charge  is  possible  with  high  quality  data 

We  proposed  a  regularization  method  for  SMC  inversion  such  that  the  TMC  is  stable  for  high  quality  data 
collected  over  a  test  stand,  and  for  measurements  taken  at  mid-height  above  a  magnetic  soil.  As  mentioned 
above,  the  SMC  method  is  data  dependent.  Therefore,  simulations  such  as  Monte  Carlo  and  inversion  of 
lesser-quality  field  data  would  be  useful  for  covering  a  wider  range  of  conditions  and  determining  the  data 
fidelity  required  for  this  method  to  succeed. 

•  Time  decay  of  Total  Magnetic  Charge  is  ordnance  specific 

Every  type  of  ordnance  and  cylinder  that  we  have  tested  so  far  shows  a  different  signature.  This  is  confirmed 
by  preliminary  canonical  analysis  and  suggests  the  high  potential  for  discrimination  of  UXO  through  their 
TMC.  In-depth  review  of  available  statistical  classification  methods  and  development  of  procedures  to  deal 
with  incomplete  recovery  of  the  TMC  curve  are  warranted. 

•  Sensitivity  to  location-orientation  parameters 

Preliminary  tests  suggest  that  the  correct  TMC  can  be  inferred  even  when  the  modeling  spheroid  has  an 
inaccurate  horizontal  position  (<  10  cm)  or  with  large  errors  in  azimuthal  direction.  Depth  error  appears 
to  affect  the  TMC  in  a  predictable  and  consistent  manner  at  all  time  channels,  thus  helping  with  the  identi¬ 
fication  of  buried  ordnance  despite  modeling  positional  error.  Detailed  sensitivity  analysis  is  warranted  to 
systematically  explore  this  effect  and  that  of  data  positional  error. 

•  If  target  depth  is  known,  a  sphere  or  a  plane  can  be  used  to  determine  the  total  magnetic  charge 

There  is  a  tremendous  amount  of  flexibility  in  the  SMC  model.  As  such,  we  have  found  that  the  shape  of  the 
magnetic  charge  surface  does  not  need  to  conform  closely  to  the  target,  i.e.,  multiple  spheroids  smaller  or 
larger  than  the  target  can  be  used  to  model  charges.  Results  not  presented  here  also  show  that  if  we  choose 
a  sphere  or  a  simple  horizontal  plane  to  model  the  charge  distribution  and  determine  the  TMC,  we  can  (1) 
still  fit  the  data  well,  and  (2)  eliminate  the  need  to  determine  the  target  orientation. 

•  Inversion  Algorithms  for  Location  and  orientation  needs  further  development 

In  this  study  the  location  and  orientation  of  the  target  was  given,  although  simple  algorithms  were  tested  to 
optimize  the  horizontal  position  to  correct  for  the  discrepancy  between  the  center  of  a  ordnance  as  defined 
in  the  field  and  the  center  of  the  modeling  spheroid  that  minimizes  data  misfit.  A  procedure  to  deal  with 
the  non-linear  part  of  data  inversion,  i.e.,  identification  of  position  and  orientation,  remains  to  be  developed 
and  tested.  At  this  point  it  also  remains  unclear  whether  data  are  needed  at  two  elevation  to  accurately 
determine  all  UXO  parameters,  as  in  Shubitidze  et  al.  (2005),  or  whether  full  non-linear  inversion  of  data 
at  a  single  elevation  would  be  possible. 


4.4  Outlook 


The  work  presented  here  has  demonstrated  that  both  the  SEA/RSS  and  SMC  approaches  have  potential  to  improve 
UXO  discrimination.  However,  a  substantial  research  and  development  effort  must  be  undertaken  before  these 
techniques  can  be  applied  to  the  real  world  UXO  problem.  We  presented  some  of  the  concerns  specific  to  the 
performance  and  suitability  of  each  method.  Addressing  these  concerns,  the  next  logical  steps  would  include: 

•  Comparison  to  Dipole  Models 

The  dipole  model  is  the  most  commonly  used  technique  for  UXO  discrimination  problems.  In  this  study, 
we  did  not  investigate  if  the  SMC  and  SEA  methods  would  represent  a  significant  improvement  to  dipole 


based  methods.  The  SMC  and  SEA  can  model  the  non-dipolar  components  of  the  secondary  field.  This 
suggests  that  the  spread  of  feature  parameters  for  a  SMC/SEA  model  would  be  relatively  smaller  than  for 
dipole  parameters.  We  would  expect  that  for  a  library  based  discrimination  approach  the  RSS  would  be 
more  suitable  than  the  dipole  model  due  to  the  ability  to  more  accurately  model  the  data.  The  SMC  model, 
however,  is  very  flexible  and,  since  a  significant  amount  of  regularization  needs  to  be  injected  into  the 
problem,  the  inversion  for  a  charge  distribution  is  much  less  straightforward  than  a  dipole  inversion.  Now 
that  a  method  for  recovering  the  charge  distribution  has  been  established,  comparisons  to  the  dipole  models 
discrimination  ability  can  be  carried  out. 

•  Sensitivity  analysis 

For  both  the  SEA  and  SMC  techniques  it  is  critical  to  determine  the  data  fidelity  required  for  application 
to  practical  discrimination  purposes.  For  the  SEA,  we  need  to  determine  a  threshold  for  the  misfit  (or 
correlation)  to  the  measured  data  below  (or  above)  which  we  classify  anomalies  as  not  coming  from  targets 
within  our  library.  For  the  SMC,  analysis  has  to  be  completed  that  determines  the  quality  of  the  total 
magnetic  charge  estimates  as  a  function  of  signal  to  noise  ratio  as  well  as  data  coverage,  positional  error, 
and  other  survey/instrument  characteristics. 

•  Multiple  Targets 

The  data  considered  in  this  study  were  from  single,  isolated  targets.  Provided  that  the  response  of  two  UXO 
targets  are  linear,  the  RSS  can  easily  model  the  total  response,  and  interpretation  techniques  based  on  the 
RSS  can  be  developed.  Nevertheless,  the  overlapping  scenario  of  UXO  plus  clutter  would  be  challenging 
to  the  SEA  techniques.  The  ability  for  the  SMC  to  be  effective  in  a  multiple  target  scenario  is  an  open 
question. 

•  Geologic  background 

A  geologic  background  signal  prevented  us  from  processing  data  collected  at  the  Sky  Research  UXO  Test 
Site.  Since  the  data  were  collected  in  a  cued  interrogation  style  (i.e.,  a  small  2  m  x  2  m  area),  we  need  to 
incorporate  a  soil  background  signal  into  our  analysis. 

We  believe  that  if  the  above  research  can  be  carried  out,  then  the  full  potential  of  the  SEA  and  SMC  tech¬ 
niques  can  be  realized  in  UXO  discrimination  surveys. 


89 


Bibliography 


Anderson,  W.  L.:  1982,  Fast  hankel  transforms  using  related  and  lagged  convolutions,  ACM  Trans.  Math.  Softw. 
8(4),  344-368. 

Andreasen,  M.  A.:  1965,  Scattering  from  bodies  of  revolution,  IEEE  Transactions  on  Antennas  and  Propagation 
13(2),  303-310. 

Ao,  C.  O.,  Braunisch,  H.,  O’Neill,  K.  and  Kong,  J.:  2002,  Quasi-magnetostatic  solution  for  a  conducting  and 
permeable  spheroid  object  to  arbitrary  excitation,  IEEE  Transactions  Geoscience  &  Remote  Sensing  40,  887- 
897. 

Arfken,  G.  B.:  1995,  Mathematical  Methods  for  Physicists,  fourth  edn.  Academic  Press,  San  Diego,  CA. 

Barrows,  B.  E.,  O’Neill,  K.,  Grzegorczyk,  T.,  Chen,  X.  and  Kong,  J.  A.:  2004,  Broadband  analytical  magneto- 
quasi-static  electromagnetic  induction  solution  for  a  conducting  and  permeable  spheroid,  IEEE  Transactions 
on  Geoscience  and  Remote  Sensing  42(11),  2479-2489. 

Beran,  L.:  2005,  Classification  of  unexploded  ordnance.  Master’s  thesis.  University  of  British  Columbia. 

Braunisch,  H.,  Ao,  C.,  O’Neil,  K.  and  Kong,  J.:  2001,  Magneto-quasistatic  response  of  conducting  and  permeable 
prolate  spheroid  under  axial  excitation,  IEEE  Transactions  on  Geoscience  and  Remote  Sensing  39,  1352-1361. 

Carin,  L.,  Yu,  H.,  Dalichaouch,  Y.,  Perry,  A.  R.,  Czipott,  P.  V.  and  Baum,  C.  E.:  2001,  On  the  wideband  emi 
response  of  a  rotationally  symmetric  permeable  and  conducting  target,  IEEE  Transactions  on  Geoscience  and 
Remote  Sensing  39(6),  1206-1213. 

Chen,  X.,  O’Neill,  K.,  Barrowes,  B.,  Grzegorczyk,  T.  and  Kong,  J.:  2004,  Application  of  a  spheroidal  mode 
approach  with  differential  evolution  in  inversion  of  magneto-quasistatic  data  for  uxo  discrimination.  Inverse 
Problems  20(6),  27^-0. 

Farquharson,  C.  G.  and  Oldenburg,  D.  W.:  2004,  A  comparison  of  automatic  techniques  for  estimating  the  regu¬ 
larization  parameter  in  nonlinear  inverse  problems.  Geophysical  Journal  International  156,  41 1 — 425. 

Gao,  P,  Collins,  L.,  Garber,  P„  Geng,  N.  and  Carin,  L.:  2000,  Classification  of  landmine-like  metal  targets  using 
wideband  electromagnetic  induction,  IEEE  Transactions  Geoscience  and  Remote  Sensing  38(3),  1352-1361. 

Geng,  N„  Baum,  K.  and  Carin,  L.:  1999,  On  the  low  frequency  natural  responses  of  conducting  and  permeable 
target,  IEEE  Transactions  Geoscience  and  Remote  Sensing  37,  347-359. 

Hansen,  P.  C.:  1997,  Rank-deficient  and  discrete  ill-posed  problems:  numerical  aspects  of  linear  inversion,  SIAM, 
Philadelphia. 

Kaufman,  A.  A.  and  Keller,  G.  V.:  1985,  Inductive  Mining  Prospecting,  Elsevier. 

Menke,  W.:  1989,  Geophysical  Data  Analysis:  Discrete  Inverse  Theory,  Revised  Edition,  Academic  Press,  Inc., 
New  York. 


90 


Miller,  J.,  Bell,  T.,  Soukup,  J.  and  Keiswetter,  D.:  2001,  Simple  phenomenological  models  for  wideband 
frequency-domain  electromagnetic  induction,  IEEE  Transactions  Geoscience  and  Remote  Sensing  39(6),  1294- 
1298. 

Newman,  G.  A.,  Hohmann,  G.  W.  and  Anderson,  W.  L.:  1986,  Transient  electromagnetic  response  of  a  three- 
dimensional  body  in  a  layered  earth.  Geophysics  51(8),  1608-1627. 

Parker,  R.  L.:  1994,  Geophysical  Inverse  Theory,  Princeton  Univ.  Press. 

Pasion,  L.  R.,  Billings,  S.  D.,  Oldenburg,  D.  W.  and  Walker,  S.  E.:  2006,  Application  of  library  based  method  to 
time  domain  electromagnetic  data  for  the  identification  of  unexploded  ordnance.  Journal  of  Applied  Geophysics 
.  Accepted  for  publication. 

Pasion,  L.  R.  and  Oldenburg,  D.  W.:  2001,  A  discrimination  algorithm  for  uxo  using  time  domain  electromagnet¬ 
ics,  Journal  of  Engineering  and  Environmental  Geophysics  6(2),  91-102. 

Riggs,  L.,  Mooney,  J.  and  Lawrence,  D.:  2001,  Identification  of  metallic  mine-like  objects  using  low  frequency 
magnetic  fields,  IEEE  Transactions  on  Geoscience  and  Remote  Sensing  39(1). 

Scales,  J.  A.  and  Tenorio,  L.:  2001,  Prior  information  and  uncertainty  in  inverse  problems.  Geophysics  66(2),  389- 
397. 

Sebak,  A.,  Shafai,  L.  and  Das,  Y.:  1991,  Near-zone  fields  scattered  by  three-dimensinal  highly  conducting  perme¬ 
able  objects  in  the  field  of  an  arbitrary  loop,  IEEE  Transactions  Geoscience  and  Remote  Sensing  29,  9-15. 

Shubitidze,  F.,  ONeill,  K.,  Haider,  S.  A.,  Sun,  K.  and  Paulsen,  K.  D.:  2002,  Application  of  the  method  of  auxiliary 
sources  for  electromagnetic  induction  problem,  IEEE  Trans.  Geoscience  &  Remote  Sensing  40(4),  928-942. 

Shubitidze,  F.,  O’Neill,  K.,  Shamatava,  I.,  Sun,  K.  and  Paulsen,  K.:  2003,  Analysis  of  emi  scattering  to  sup¬ 
port  uxo  discrimination:  Heterogeneous  and  multi  objects,  SPIE  Conference  on  Detection  and  Remediation 
Technologies  for  Mines  and  Minelike  Targets,  The  International  Society  for  Optical  Engineering,  pp.  928-939. 

Shubitidze,  F.,  O’Neill,  K.,  Shamatava,  I.,  Sun,  K.  and  Paulsen,  K.:  2004,  Use  of  standardized  source  sets  for 
enhanced  emi  classification  of  buried  heterogeneous  objects,  SPIE  Conference  on  Detection  and  Remediation 
Technologies  for  Mines  and  Minelike  Targets,  The  International  Society  for  Optical  Engineering,  pp.  263-274. 

Shubitidze,  F.,  ONeill,  K.,  Shamatava,  I.,  Sun,  K.  and  Paulsen,  K.:  2005a,  A  simple  magnetic  charge  model  for 
classification  of  multiple  buried  metallic  objects  in  cases  with  overlapping  signals.  Proceedings  from  SAGEEP 
03. 

Shubitidze,  F.,  ONeill,  K.,  Shamatava,  I.,  Sun,  K.  and  Paulsen,  K.  D.:  2005b,  Fast  and  accurate  calculation 
of  physically  complete  emi  response  by  a  heterogeneous  metallic  object,  IEEE  Trans.  Geoscience  &  Remote 
Sensing  43(8),  1736-1750. 

Shubitidze,  F.,  O’Neill,  K.,  Sun,  K.  and  Paulsen,  K.:  2004,  Investigation  of  broadband  electromagnetic  induction 
scattering  by  highly  conducting,  permeable,  arbitrarily  shaped  3-d  objects,  IEEE  Transactions  Geoscience  and 
Remote  Sensing  42(3),  928-942. 

Shubitidze,  F.,  O’Neill,  K.,  Sun,  K.  and  Shamatava,  I.:  2004,  Coupling  between  highly  conducting  and  permeable 
metallic  objects  in  the  emi  frequency  range.  Applied  Computational  Electromagnetic  Society  Journal,  Special 
issue  on  ACES  2003  conference,  Vol.  19,  pp.  139-148. 

Shubitidze,  F.,  O’Neill,  K.,  Sun,  K.,  Shamatava,  I.  and  Paulsen,  K.:  2004,  A  hybrid  full  MAS  and  combined 
MAS/TS A  algorithm  for  electromagnetic  induction  problems  sensing.  Applied  Computational  Electromagnetic 
Society  Journal,  Special  issue  on  ACES  2003  conference,  Vol.  19,  pp.  112-127. 

Sun,  K.,  O’Neill,  K.,  Shubitidze,  F.,  Haider,  S.  A.  and  Paulsen.,  K.  D.:  2002,  Simulation  of  electromagnetic 
induction  scattering  from  targets  with  negligible  to  moderate  penetration  by  primary  fields,  IEEE  transactions 
on  geoscience  and  remote  sensing  40,  910-927. 


91 


Sun,  K.,  O’Neill,  K.,  Shubitidze,  F.,  Shamatava,  I.  and  Paulsen.,  K.  D.:  2004,  Theoretical  analyses  of  TSA 
formulation  and  its  domain  of  validity,  IEEE  transactions  on  geoscience  and  remote  sensing  42(9),  1871-1881. 

Sun,  K.,  O’Neill,  K.,  Shubitidze,  F.,  Shamatava,  I.  and  Paulsen.,  K.  D.:  2005,  Fast  data-derived  fundamental 
spheroidal  excitation  models  with  application  to  uxo  identification,  IEEE  transactions  on  geoscience  and  re¬ 
mote  sensing  43(11),  2573-2583. 

Tarantola,  A.:  1987,  Inverse  Problem  Theory:  Methods  for  Data  Fitting  and  Model  Parameter  Estimation,  Else¬ 
vier. 

Zhang,  Y.,  Collins,  L.,  Yu,  H.,  Baum,  C.  E.  and  Carin,  L.:  2003,  Sensing  of  unexploded  ordnance  with  mag¬ 
netometer  and  induction  data:  Theory  and  signal  processing,  IEEE  Transactions  on  Geoscience  and  Remote 
Sensing  41(5),  1005-1015. 


92 


